Department Behavioural Neurobiology, Max Planck Institute for Biological Intelligence
Coucal Project
Department Biology II, Ludwig Maximilians University Munich
Published
November 24, 2023
In this document we provide all the necessary code for reproducing the analyses presented in our manuscript. To access the dataset and Rmarkdown file, please download this GitHub repository. An explanation of the files in the repository can be found in the Readme file. Please don’t hesitate to contact Luke at (luke.eberhart[at]bi.mpg.de) if you have any questions regarding the code below. For questions regarding the data collection and the study system, please contact Wolfgang (wgoymann[at]bi.mpg.de) or Safari (safariignas[at]yahoo.co.uk)
Prerequisites
GitHub repository
Simply follow the link to the GitHub and click on Download ZIP on the right-hand side of the page. After unzipping the folder, open the coucal_demography.Rproj file in RStudio, which will initiate a new RStudio session with the project directory referencing all the other files downloaded in the zipped folder (e.g., output/bootstraps/sensitivity_analysis/.
R packages
The following packages are needed for analysis and can be easily installed from CRAN or GitHub by running the following code chunk:
Code
# a vector of all the packages needed in the projectpackages_required_in_project <-c("RMark","tidyverse","readxl","BaSTA","pbapply","RColorBrewer","grid","Rmisc","gss","arm","partR2","parameters","standardize","colorBlindness","ggthemes","patchwork","gt","rptR","tidybayes","broom.mixed","effects","patchwork","devtools","unmarked","R2ucare","marked","merTools","bootpredictlme4","extrafont","survminer","bdsmatrix","coxme","epitools","survival","magrittr","ggpubr","reshape2","sp","adehabitatLT","reshape","coefplot","mapview","lubridate","effects","merDeriv","remotes")# of the required packages, check if some need to be installednew.packages <- packages_required_in_project[!(packages_required_in_project %in%installed.packages()[,"Package"])]# install all packages that are not locally availableif(length(new.packages)) install.packages(new.packages)# install the bootpredictlme4 package directly from GitHubremotes::install_github("RemkoDuursma/bootpredictlme4")# load all the packages into the current R sessionlapply(packages_required_in_project, require, character.only =TRUE)
Plotting themes
The following plotting themes, colors, and typefaces are used throughout the project:
Code
# Find fonts from computer that you want. Use regular expressions to do this# For example, load all fonts that are 'candara' or 'Candara'extrafont::font_import(pattern ="[V/v]erdana", prompt =FALSE) # check which fonts were loadedextrafont::fonts()extrafont::fonttable()extrafont::loadfonts() # load these into R# define the plotting theme to be used in subsequent ggplotsluke_theme <-theme_bw() +theme(text =element_text(family ="Verdana"),legend.title =element_text(size =16),legend.text =element_text(size =12),axis.title.x =element_text(size =12),axis.text.x =element_text(size =10), axis.title.y =element_text(size =12),axis.text.y =element_text(size =10),strip.text =element_text(size =12),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_line(size =0.5, colour ="grey40"),axis.ticks.length =unit(0.2, "cm"),panel.border =element_rect(linetype ="solid", colour ="grey"),legend.position =c(0.1, 0.9) )# set plotting color palettessex_pal2 <-c(pull(ggthemes_data$wsj$palettes$colors6[3,2]), pull(ggthemes_data$wsj$palettes$colors6[2,2]))sex_pal3 <-c(pull(ggthemes_data$wsj$palettes$colors6[3,2]), pull(ggthemes_data$wsj$palettes$colors6[3,2]),pull(ggthemes_data$wsj$palettes$colors6[2,2]),pull(ggthemes_data$wsj$palettes$colors6[2,2]))# specify the facet labels for each species and sexspecies_names <-c('BC'="Black coucal",'WBC'="White-browed coucal")sex_names <-c('female'="Females",'male'="Males")analysis_names <-c('male'="Male Mo scenario",'female'="Female Mo scenario")# color of mean estimate point in forest plotscol_all <-"#2E3440"# custom color palette for the plotting of Juvenile and Adult statscbPalette_LTRE <-c("#D9D9D9", "#D9D9D9", "#D9D9D9", "#D9D9D9", "#A6A6A6", "#A6A6A6","#A6A6A6")cbPalette_sex_diff <-c("#D9D9D9", "#D9D9D9", "#D9D9D9", "#D9D9D9", "#A6A6A6")# plot the comparative LTRE resultsvital_rate_theme <-theme_bw() +theme(text =element_text(family ="Verdana"),legend.position ="none",panel.grid.major =element_blank(),panel.grid.minor =element_blank(),axis.ticks.length =unit(0.1, "cm"),panel.border =element_blank(),panel.spacing.x =unit(0.3, "lines"),panel.spacing.y =unit(0.7, "lines"),strip.background =element_blank() )species.labs <-c("Black Coucal", "White-browed Coucal")names(species.labs) <-c("BC", "WBC")
Analytical Functions
The following custom functions are used to bootstrap the various datasets and run the stochastic demographic model. To reproduce the results presented in our manuscript, these functions need to be stored in the environment of your R Session.
adult_CJS_bootstrap_data()
This function randomly samples the adult mark-recapture dataset to produce a random one the same size as the original, with replacement
arguments:
adult: the adult mark-recapture dataset
num_boot: a unique number for the pbapply simulation (don’t specify)
species: species name (either “BC” or “WBC”)
iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
Code
adult_CJS_bootstrap_data <-function(adult, num_boot, species, iter_add) {# sample a new adult mark-recapture dataset the same size as the original, # with replacement adult_boot <- adult %>%sample_n(nrow(.), replace =TRUE)# make a list of the new dataset, iteration, species out <-list(adult_boot = adult_boot,iter = num_boot + ((iter_add -1) * niter),species = species)}
bootstrap_adult_survival_CJS()
runs the CJS survival analyses using the the bootstrapped sample created from adult_CJS_bootstrap_data().
arguments:
coucal_boot_list: the bootstrapped adult mark-recapture dataset produced from the adult_CJS_bootstrap_data() function above
num_boot: a unique number for the pbapply simulation (don’t specify)
first_year: first year of the study (i.e., the year of the first event in the encounter history)
bootstrap_name: a unique string for naming of files in the output directory
species: species name (either “BC” or “WBC”)
iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
Code
bootstrap_adult_survival_CJS <-function(coucal_boot_list, num_boot, first_year, bootstrap_name, species, iter_add, prefix_number,out_directory ="/output/bootstraps/raw/") {# specify the bootstrapped data samples (from the previous function) adult_ch <- coucal_boot_list[["adult_boot"]]# clean up capture histories adult_ch <- adult_ch %>%ungroup() %>%as.data.frame() %>%mutate(sex =as.factor(sex))# process the adult data as a "CJS" analysis coucal_adult.proc <- RMark::process.data(adult_ch, model ="CJS",groups =c("sex"),begin.time = first_year)# Create the design matrix from the processed mark-recapture datasets coucal_adult.ddl <- RMark::make.design.data(coucal_adult.proc) adult_survival_analysis_run =function(proc_data, design_data){# sex- and stage-specific survival: Phi.sex =list(formula =~ sex) # Models exploring variation in encounter probability# constant: p.dot =list(formula =~1)# sex-dependent: p.sex =list(formula =~ sex)# factorial variation across year: p.Time =list(formula =~ time)# additive effects of sex and factorial year: p.sex_Time =list(formula =~ sex + time)# create a list of candidate models for all the a models above that begin with # either "Phi." or "p." cml <- RMark::create.model.list("CJS")# specify the data, design matrix, delete unneeded output files, and # run the models in Program MARK model.list <- RMark::mark.wrapper(cml, data = proc_data, ddl = design_data, delete =TRUE, wrap =FALSE, threads =1, brief =TRUE,silent =TRUE, output =FALSE, prefix = prefix_number)# output the model list and sotre the resultsreturn(model.list) }# Run the models on the bootstrapped data adult_survival_analysis_out <-adult_survival_analysis_run(proc_data = coucal_adult.proc,design_data = coucal_adult.ddl) extract_top_model_output <-function(rmark_output, top_model =TRUE, mod_num){# Find the model number for the first ranked model of the AIC tableif(top_model ==TRUE){ mod_num <-as.numeric(rownames(rmark_output$model.table[1,])) }else{ mod_num <- mod_num }# extract and wrangle reals from model output reals <- rmark_output[[mod_num]]$results$real %>%bind_cols(data.frame(str_split_fixed(rownames(.), " ", n =5)), .) %>%filter(X1 =="Phi") %>%mutate(sex =as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) =="F", "Female", "Male"))) %>% dplyr::select(sex, estimate) %>%mutate(iter = num_boot + ((iter_add -1) * niter),species = species)# extract and wrangle the beta estimates from the model output betas <- rmark_output[[mod_num]]$results$beta %>%mutate(statistic =rownames(.),iter = num_boot + ((iter_add -1) * niter),species = species) %>%mutate(parameter =ifelse(grepl(x = statistic, pattern ="S:"), "S",ifelse(grepl(x = statistic, pattern ="p:"), "p",ifelse(grepl(x = statistic, pattern ="r:"), "r",ifelse(grepl(x = statistic, pattern ="F:"), "F", ifelse(grepl(x = statistic, pattern ="Phi:"), "Phi", "XXX"))))),variable =ifelse(grepl(x = statistic, pattern ="Intercept"), "Intercept",ifelse(grepl(x = statistic, pattern ="sexM:Cubic"), "sexM:Cubic",ifelse(grepl(x = statistic, pattern ="sexM:Quadratic"), "sexM:Quadratic",ifelse(grepl(x = statistic, pattern ="sexM:Time"), "sexM:Time",ifelse(grepl(x = statistic, pattern ="sexM"), "sexM",ifelse(grepl(x = statistic, pattern ="Time"), "Time",ifelse(grepl(x = statistic, pattern ="Cubic"), "Cubic", ifelse(grepl(x = statistic, pattern ="Quadratic"), "Quadratic","XXX"))))))))) %>% dplyr::select(-statistic)# consolidate the AIC model selection results AIC_table <- rmark_output$model.table %>%mutate(iter = num_boot + ((iter_add -1) * niter),species = species,model_no_orig =as.numeric(rownames(.))) %>%mutate(model_no_rank =as.numeric(rownames(.)))# consolidate the all output into a list survival_model_output_list <-list(reals = reals,betas = betas,AIC_table = AIC_table) survival_model_output_list }# extract and format survival rates from model output adult_survival_model_output_list <-extract_top_model_output(rmark_output = adult_survival_analysis_out)# make a list of all the results from this iteration bootstrap_results_list <-list(adult_survival_model_output = adult_survival_model_output_list,bootstrapped_data = coucal_boot_list)# extract file directory file_directory <-paste0(getwd(), "/", out_directory, "/", bootstrap_name, "_", num_boot, ".Rds")# save the results of the iterationsave(bootstrap_results_list, file = file_directory) bootstrap_results_list }
run_bootstrap_adult_survival_CJS()
run the bootstrap and the CJS analysis (passed through the pbapply::pbsapply() function), see arguments defined above
arguments:
adult: the adult mark-recapture dataset
num_boot: a unique number for the pbapply() simulation (don’t specify)
species: species name (either “BC” or “WBC”)
iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
first_year: first year of the study (i.e., the year of the first event in the encounter history)
bootstrap_name: a unique string for naming of files in the output directory
prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
Code
run_bootstrap_adult_survival_CJS <-function(adult, num_boot, species, iter_add, first_year, bootstrap_name, prefix_number, out_directory){# run the sampling function and specify the datasets bootstrap_data_list <-adult_CJS_bootstrap_data(adult = adult, num_boot = num_boot,species = species, iter_add = iter_add)# run the survival analysis and ASR deduction on the sampled data result <-bootstrap_adult_survival_CJS(coucal_boot_list = bootstrap_data_list, num_boot = num_boot, first_year = first_year,bootstrap_name = bootstrap_name,species = species,iter_add = iter_add,prefix_number = prefix_number,out_directory = out_directory) result}
surv_boot_out_wrangle()
loads all the output rds files from the parallel adult CJS bootstrapping procedure and consolidates them into a single object for further analysis
arguments:
species: species name (either “BC” or “WBC”)
niter: number of iterations specified in the stochastic simulation (used to predict how many files to merge from the output directory)
out_directory: location in the project directory where the output files from the simulation were stored (e.g., “tmp”)
Code
# boot_out_wrangle() surv_boot_out_wrangle <-function(species, niter, output_directory ="output/bootstraps/"){# specify directory string boot_out_path <-paste0(output_directory, species,"_adult_survival_CJS_bootstrap_result.rds")# load the bootstrap output file survival_bootstrap_result <-readRDS(boot_out_path)# bind all the adult rates output together adult_reals_survival_rates_boot <-rbind(do.call(rbind, lapply(seq(from =1, to = niter, by =1),function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$reals))) %>%arrange(iter, sex)# bind all the adult beta estimates output together adult_betas_survival_rates_boot <-rbind(do.call(rbind, lapply(seq(from =1, to = niter, by =1), function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$betas)))# bind all the adult AIC table output together adult_AIC_tables_boot <-rbind(do.call(rbind, lapply(seq(from =1, to = niter, by =1), function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$AIC_table)))# tidy all the bound data together into a mega list of results boot_out_tidy <-list(adult_reals_survival_rates_boot = adult_reals_survival_rates_boot,adult_betas_survival_rates_boot = adult_betas_survival_rates_boot,adult_AIC_tables_boot = adult_AIC_tables_boot) boot_out_tidy}
approx_sd()
derive the approximate sd given the upper and lower 95% confidence intervals (see:Cochrane Chapter 7.7.3.2)
arguments:
x1: lower bound of the confidence interval
x2: upper bound of the confidence interval
Code
# function to approximate the sd given the 95% CIsapprox_sd <-function(x1, x2){ (x2-x1) / (qnorm(0.975) -qnorm(0.025))}
matrix_ASR()
calculates the ASR of the population based on the two-sex two-stage projection matrix built by the coucal_matrix() function (below).
arguments:
M: a two-sex x by x projection matrix produced by the coucal_matrix() function above
n: x-lengthed vector representing starting stage distribution (the default is a vector with 10 individuals in each stage)
h: harem size index (default is 1, monogamy; polyandry < 1, polygyny > 1)
k: clutch size (default is 4 eggs, mean for both coucal species)
iterations: number of time steps to project the demogrphic model (n.b., after 100 steps our model reaches the stable stage distribution)
HSR: hatching sex ratio (expressed as the proportion of male chicks)
num_boot: a unique number for the pbapply simulation (don’t specify)
species: species name (either “BC” or “WBC”)
iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk (if you run this on one session, then don’t specify)
Code
matrix_ASR <-function(M, n =rep(10, nrow(M)), h =1, k =4, iterations =1000, HSR =0.5, num_boot, species, iter_add){# Number of stages in matrix x <-length(n)# Number of time steps to simulate t <- iterations# an empty t by x matrix to store the stage distributions stage <-matrix(data =numeric(x * t), nrow = x, ncol = t)# an empty t by x matrix to store the stage-specific frequencies pop_dist <-matrix(data =numeric(x * t), nrow = x, ncol = t)rownames(pop_dist) <-rownames(M)colnames(pop_dist) <-0:(t -1)# an empty t vector to store the population sizes pop <-numeric(t)# for loop that goes through each of t time stepsfor (i in1:t) {# stage distribution at time t stage[,i] <- n# population size at time t pop[i] <-sum(n)# number of male adults at time t = (number alive at t-1) * (survival rate) M2 <- stage[4, i] * M["M_Adult", "M_Adult"]# number of female adults at time t F2 <- stage[2, i] * M["F_Adult", "F_Adult"]# Female freq-dep fecundity of Female chicks M["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1- HSR))# Female freq-dep fecundity of Male chicks M["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)# Male freq-dep fecundity of Female chicks M["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1- HSR))# Male freq-dep fecundity of Male chicks M["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)# define the new n (i.e., new stage distribution at time t) n <- M %*% n# define rownames of stage matrixrownames(stage) <-rownames(M)# define colnames of stage matrixcolnames(stage) <-0:(t -1)# calculate the proportional stable stage distribution stage <-apply(stage, 2, function(x) x/sum(x))# define stable stage as the last stage stable.stage <- stage[, t] pop_dist[, i] <- n }# calc ASR as the proportion of the adult stable stage class that is male ASR <- stable.stage[4]/(stable.stage[2] + stable.stage[4]) SSR <- stable.stage[3]/(stable.stage[1] + stable.stage[3])# make a list of results pop.proj <-list(ASR = ASR,SSR = SSR,lambda = pop[t]/pop[t -1],pop_size = pop,stage_size = pop_dist,stable.stage = stable.stage,stage.vectors = stage,SSD_M2 = stable.stage[4],SSD_F2 = stable.stage[2],iter = num_boot + ((iter_add -1) * niter), species = species)# print the list as output to the function pop.proj }
coucal_matrix()
builds the two-sex Lefkovitch matrix using the vital rates specified in the demographic_rates object. This function also give the option to specify a one-sex matrix for analytical comparisons between 2- and 1-sex matrix models
arguments:
demographic_rates_: a list of all the vital rates of an iteration of the stochastic simulation run within the run_bootstrap_juv_hazd_ad_surv_ASR() function below
two_sex: whether to make a one-sex or two-sex matrix of the rates provided (default is ‘two_sex = TRUE’)
Code
coucal_matrix <-function(demographic_rates_, two_sex =TRUE){if(two_sex){# Define coucal life-stages of the coucal matrix model stages <-c("F_1st_year", "F_Adult", "M_1st_year", "M_Adult")# Build the 4x4 matrix result <-matrix(c(# top row of matrix0, NA, 0, NA, # second row of matrix (demographic_rates_$Egg_survival * demographic_rates_$F_Nestling_survival * demographic_rates_$F_Groundling_survival * demographic_rates_$F_Fledgling_survival), demographic_rates_$F_Adult_survival,0, 0,# third row of matrix0, NA, 0, NA, # fourth row of matrix0, 0, (demographic_rates_$Egg_survival * demographic_rates_$M_Nestling_survival * demographic_rates_$M_Groundling_survival * demographic_rates_$M_Fledgling_survival), demographic_rates_$M_Adult_survival),nrow =length(stages), byrow =TRUE,dimnames =list(stages, stages)) }else{# Define coucal life-stages of the Ceuta snowy coucal matrix model stages <-c("1st_year", "Adult")# Build the 4x4 matrix result <-matrix(c(# top row of matrix0, NA,# second row of matrix (demographic_rates_$Egg_survival * demographic_rates_$Nestling_survival * demographic_rates_$Groundling_survival * demographic_rates_$Fledgling_survival), demographic_rates_$Adult_survival),nrow =length(stages), byrow =TRUE,dimnames =list(stages, stages)) } result }
bootstrap_hazard_data()
This function randomly samples one sibling per nest from the offspring survival dataset
arguments:
offspring: the radio-tracked offspring dataset
num_boot: a unique number for the pbapply simulation (don’t specify)
species: species name (either “BC” or “WBC”)
iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd() function (n.b., Gu 2014 J. Stat Software 58:art5 recommends alpha = 1.4)
max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
niter: number of iterations specified in the stochastic simulation
Code
# bootstrap_hazard_data() bootstrap_hazard_data <-function(offspring, num_boot, species, iter_add, alpha_value, max_time, niter) {# set attempt to 0 at start of each loop attempt <-0 max_time <- max_time hzd_ss_try_M_lambda <-1.1 hzd_ss_try_F_lambda <-1.1# store simulated estimates only if lambda is less than 1.1while( (hzd_ss_try_M_lambda >-1.1| hzd_ss_try_F_lambda >1.1) && attempt <=100 ) {# next attempt attempt <- attempt +1# sample a new offspring dataset containing only one nest member per draw offspring_boot <- offspring %>%group_by(nest_ID) %>%sample_n(1)try( hzd_ss_try_M <-sshzd(Surv(exit, event, entry) ~ exit,data =filter(offspring_boot, sex =="M"),alpha = alpha_value),silent =TRUE )try( hzd_ss_try_F <-sshzd(Surv(exit, event, entry) ~ exit,data =filter(offspring_boot, sex =="F"),alpha = alpha_value),silent =TRUE )# store calculated peak (i.e., the apex of the polynomial curve)try( hzd_ss_try_M_lambda <- hzd_ss_try_M$lambda )# store calculated peak (i.e., the apex of the polynomial curve)try( hzd_ss_try_F_lambda <- hzd_ss_try_F$lambda ) }# store simulated estimates only if the hzdcurve can be estimatedwhile( (exists("hzd_curve_try_M") ==FALSE|exists("hzd_curve_try_F") ==FALSE) && attempt <= max_time) { time_vector <-seq(0, max_time - attempt, 1)# next attempt attempt <- attempt +1try( hzd_ss_try_M <-sshzd(Surv(exit, event, entry) ~ exit, data =filter(offspring_boot, sex =="M"), alpha = alpha_value),silent =TRUE )try( hzd_ss_try_F <-sshzd(Surv(exit, event, entry) ~ exit, data =filter(offspring_boot, sex =="F"), alpha = alpha_value),silent =TRUE )# simulate an estimatetry( hzd_curve_try_M <-hzdcurve.sshzd(object = hzd_ss_try_M, time = time_vector, se =TRUE),silent =TRUE )try( hzd_curve_try_F <-hzdcurve.sshzd(object = hzd_ss_try_F, time = time_vector, se =TRUE),silent =TRUE ) }# make a list including elements that will be used in the next function out <-list(offspring_boot = offspring_boot, iter = num_boot + ((iter_add -1) * niter),species = species,time_vector = time_vector,lambdaM = hzd_ss_try_M_lambda,lambdaF = hzd_ss_try_F_lambda) }
run_bootstrap_juv_hazd_ad_surv_ASR()
initiates the bootstrap_hazard_data() and runs the juvenile survival analysis to derive sex- and stage-specific survival rates. Then constructs the matrix model with the rates and projects this into the future to acquire the stable stage distribution and derive the ASR.
arguments:
num_boot: a unique number for the pbapply simulation (don’t specify)
offspring: the radio-tracked offspring dataset
k_dist: a two-number vector containing the mean and sd of clutch size (output from the analyses below)
HSR_dist: a two-number vector containing the mean and sd of hatching sex ratio (output from the analyses below)
h_dist: a two-number vector containing the mean and sd of the harem size index (output from the analyses below)
egg_surv_dist: a two-number vector containing the mean and sd of egg survival (output from the analyses below)
flight_age_distF: a two-number vector containing the mean and sd of female-specific flight age (output from the analyses below)
flight_age_distM: a two-number vector containing the mean and sd of male-specific flight age (output from the analyses below)
fledge_age_distF: a two-number vector containing the mean and sd of female-specific fledge age (output from the analyses below)
fledge_age_distM: a two-number vector containing the mean and sd of male-specific fledge age (output from the analyses below)
bootstrap_name: a unique string for naming of files in the output directory
adult_surival_boot_out: the bootstrapped adult encounter history data
species: species name (either “BC” or “WBC”)
iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd function (default is 1.4).
max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
niter: number of iterations specified in the stochastic simulation
out.directory: location in the project directory where the output files should be stored (e.g., “tmp”)
Code
run_bootstrap_juv_hazd_ad_surv_ASR <-function(num_boot, offspring, k_dist, HSR_dist, h_dist, egg_surv_dist, flight_age_distF, flight_age_distM, fledge_age_distF, fledge_age_distM, bootstrap_name, adult_surival_boot_out, species, iter_add, prefix_number,alpha_value =1.4,max_time =70, niter, output.directory){# run the sampling function and specify the datasets coucal_boot_list <-bootstrap_hazard_data(offspring = offspring, num_boot = num_boot,species = species, iter_add = iter_add,alpha_value = alpha_value,max_time = max_time,niter = niter)# run the survival analysis and ASR deduction on the sampled data# specify the bootstrapped data samples (from the previous function) offspring_data <- coucal_boot_list[["offspring_boot"]]# clean up capture histories offspring_data <- offspring_data %>%ungroup() %>%as.data.frame()# fit smoothed spline of hazard function for either sex M_haz_ss <-sshzd(Surv(exit, event, entry) ~ exit, data =filter(offspring_data, sex =="M"), alpha = alpha_value) F_haz_ss <-sshzd(Surv(exit, event, entry) ~ exit, data =filter(offspring_data, sex =="F"), alpha = alpha_value) haz_ss_lambda <-list(Male_haz_ss_lambda = M_haz_ss$lambda,Female_haz_ss_lambda = F_haz_ss$lambda)# extract fitted estimates from the spline function M_haz_ss_curve <-hzdcurve.sshzd(object = M_haz_ss, time = coucal_boot_list[["time_vector"]], se =TRUE) F_haz_ss_curve <-hzdcurve.sshzd(object = F_haz_ss, time = coucal_boot_list[["time_vector"]], se =TRUE) haz_ss_curve <-expand.grid(species = species, age = coucal_boot_list[["time_vector"]],sex =c("Male", "Female")) %>%mutate(fit =c(M_haz_ss_curve$fit, F_haz_ss_curve$fit),se =c(M_haz_ss_curve$se, F_haz_ss_curve$se)) %>%mutate(estimate =1- fit,upper =1- fit *exp(1.96* se),lower =1- fit /exp(1.96* se),iter = num_boot)# transform the daily nestling survival (DCS) to apparent fledgling success# by calculating the product of all DCS estimates: fledge_ageM <-rnorm(1, fledge_age_distM[1], fledge_age_distM[2]) fledge_ageF <-rnorm(1, fledge_age_distF[1], fledge_age_distF[2]) flight_ageM <-rnorm(1, flight_age_distM[1], flight_age_distM[2]) flight_ageF <-rnorm(1, flight_age_distF[1], flight_age_distF[2]) coucal_nestling_survivalF <- haz_ss_curve %>%filter(age <= fledge_ageF & sex =="Female") %>%group_by(sex) %>% dplyr::summarise(value =prod(estimate), .groups ='drop') %>%mutate(stage ="nestling",rate ="survival") coucal_nestling_survivalM <- haz_ss_curve %>%filter(age <= fledge_ageM & sex =="Male") %>%group_by(sex) %>% dplyr::summarise(value =prod(estimate), .groups ='drop') %>%mutate(stage ="nestling",rate ="survival") coucal_nestling_survival <-bind_rows(coucal_nestling_survivalF, coucal_nestling_survivalM) coucal_groundling_survivalF <- haz_ss_curve %>%filter(age < (flight_ageF + fledge_ageF) & age > fledge_ageF & sex =="Female") %>%group_by(sex) %>% dplyr::summarise(value =prod(estimate), .groups ='drop') %>%mutate(stage ="groundling",rate ="survival") coucal_groundling_survivalM <- haz_ss_curve %>%filter(age < (flight_ageM + fledge_ageM) & age > fledge_ageM & sex =="Male") %>%group_by(sex) %>% dplyr::summarise(value =prod(estimate), .groups ='drop') %>%mutate(stage ="groundling",rate ="survival") coucal_groundling_survival <-bind_rows(coucal_groundling_survivalF, coucal_groundling_survivalM) coucal_fledgling_survivalF <- haz_ss_curve %>%filter(age >= (flight_ageF + fledge_ageF) & sex =="Female") %>%group_by(sex) %>% dplyr::summarise(value =prod(estimate), .groups ='drop') %>%mutate(stage ="fledgling",rate ="survival") coucal_fledgling_survivalM <- haz_ss_curve %>%filter(age >= (flight_ageM + fledge_ageM) & sex =="Male") %>%group_by(sex) %>% dplyr::summarise(value =prod(estimate), .groups ='drop') %>%mutate(stage ="fledgling",rate ="survival") coucal_fledgling_survival <-bind_rows(coucal_fledgling_survivalF, coucal_fledgling_survivalM) coucal_adult_survival <-filter(adult_surival_boot_out, iter ==sample(1:niter, 1)) %>% dplyr::select(-species, -iter) %>%mutate(stage ="adult",rate ="survival") egg_survival <-rnorm(1, egg_surv_dist[1], egg_surv_dist[2]) coucal_egg_survival <-data.frame(sex =NA,value = egg_survival,stage =c("egg"),rate =c("survival")) coucal_mating_system <-data.frame(sex =NA,value =1/rnorm(1, h_dist[1], h_dist[2]),stage =c("h"),rate =c("fecundity")) coucal_clutch_size <-data.frame(sex =NA,value =rnorm(1, k_dist[1], k_dist[2]),stage =c("k"),rate =c("fecundity")) coucal_HSR <-data.frame(sex =NA,value =rnorm(1, HSR_dist[1], HSR_dist[2]),stage =c("HSR"),rate =c("fecundity")) coucal_flight_age <-data.frame(sex =c("Female", "Male"),value =c(flight_ageF, flight_ageM),stage =c("flight_age"),rate =c("development")) coucal_fledge_age <-data.frame(sex =c("Female", "Male"),value =c(fledge_ageF, fledge_ageM),stage =c("fledge_age"),rate =c("development"))# Bind the juvenile and adult dataframe with the nestlings coucal_vital_rates <-bind_rows(coucal_egg_survival, coucal_nestling_survival, coucal_groundling_survival, coucal_fledgling_survival, coucal_adult_survival, coucal_mating_system, coucal_clutch_size, coucal_HSR, coucal_flight_age, coucal_fledge_age) %>%as.data.frame() %>%mutate(iter = num_boot, #+ ((iter_add - 1) * niter),species = species) %>%arrange(sex, rate, stage)# Create a list of demographic rates from the survival analyses above demographic_rates <-list(Egg_survival =pull(filter(coucal_vital_rates, stage =="egg"), value),F_Nestling_survival =pull(filter(coucal_vital_rates, stage =="nestling"& rate =="survival"& sex =="Female"), value),F_Groundling_survival =pull(filter(coucal_vital_rates, stage =="groundling"& rate =="survival"& sex =="Female"), value),F_Fledgling_survival =pull(filter(coucal_vital_rates, stage =="fledgling"& rate =="survival"& sex =="Female"), value),F_Adult_survival =pull(filter(coucal_vital_rates, stage =="adult"& rate =="survival"& sex =="Female"), value),M_Nestling_survival =pull(filter(coucal_vital_rates, stage =="nestling"& rate =="survival"& sex =="Male"), value),M_Groundling_survival =pull(filter(coucal_vital_rates, stage =="groundling"& rate =="survival"& sex =="Male"), value),M_Fledgling_survival =pull(filter(coucal_vital_rates, stage =="fledgling"& rate =="survival"& sex =="Male"), value),M_Adult_survival =pull(filter(coucal_vital_rates, stage =="adult"& rate =="survival"& sex =="Male"), value),# Define hatching sex ratioHSR =pull(filter(coucal_vital_rates, stage =="HSR"), value),# Define the mating system (h), and clutch size (k)h =pull(filter(coucal_vital_rates, stage =="h"), value),k =pull(filter(coucal_vital_rates, stage =="k"), value))# Build matrix based on rates specified in the list above demographic_matrix <-coucal_matrix(demographic_rates_ = demographic_rates)# Determine the ASR at the stable stage distribution ASR_SSD <-matrix_ASR(M = demographic_matrix,h = demographic_rates$h,HSR = demographic_rates$HSR, iterations =100,num_boot = num_boot,species = species,iter_add =1,n =rep(10, nrow(demographic_matrix)))# Extract ASR ASR_estimate <- ASR_SSD$ASR# make a list of all the results from this iteration bootstrap_results_list <-list(offspring_hazard_lambda = haz_ss_lambda, offspring_hazard_rates = haz_ss_curve,coucal_vital_rates = coucal_vital_rates,ASR_SSD = ASR_SSD,bootstrapped_data = coucal_boot_list)# extract file directory file_directory <-paste0(getwd(), output.directory, "/", bootstrap_name, "_", num_boot, ".Rds")# save the results of the iterationsave(bootstrap_results_list, file = file_directory) result = bootstrap_results_list result}
hazard_boot_out_wrangle()
If the product of run_bootstrap_juv_hazd_ad_surv_ASR() is saved on disk, this function will tidy up the components of the raw output for subsequent analyses
arguments:
species: species name (either “BC” or “WBC”)
niter: number of iterations specified in the stochastic simulation
output_dir: location in the project directory where the output files are to be found by the function (e.g., “tmp”)
rds_file: name of the stochastic simulation rds file saved to disk
Code
# single_model_boot_out_wrangle() hazard_boot_out_wrangle <-function(species, niter, output_dir, rds_file){# specify directory string boot_path <-paste0(output_dir, species, rds_file, ".rds")# load each of the four bootstrap output files bootstrap_out <-readRDS(file = boot_path)# bind all the fledgling daily survival rates output together hazard_rates_boot <-do.call(rbind, lapply(seq(from =1, to = niter, by =1),function(x) bootstrap_out["offspring_hazard_rates", x]$offspring_hazard_rates)) %>%arrange(iter, sex, age)# bind all the ASR estimates output together ASR_ests_boot <-do.call(rbind, lapply(seq(from =1, to = niter, by =1),function(x) bootstrap_out["ASR_SSD", x]$ASR_SSD$ASR)) %>%as.data.frame() %>%mutate(species = species)# bind all the ASR estimates output together vital_rate_ests_boot <-do.call(rbind, lapply(seq(from =1, to = niter, by =1),function(x) bootstrap_out["coucal_vital_rates", x]$coucal_vital_rates))# tidy all the bound data together into a mega list of results boot_out_tidy <-list(hazard_rates_boot = hazard_rates_boot,vital_rate_ests_boot = vital_rate_ests_boot,ASR_ests_boot = ASR_ests_boot) boot_out_tidy }
sex_diff_hazard()
takes the list of results from the bootstrap procedure and calculates the sex difference in survival at each life stage for each iteration
arguments:
boot_out_list: the object produced by the hazard_boot_out_wrangle() above
niter: number of iterations specified in the stochastic simulation
Code
sex_diff_hazard <-function(boot_out_list, niter) {# make an empty datarame to store the results sex_diff_surv_output <-data.frame(Adult = niter,Fledgling = niter,Groundling = niter,Nestling = niter,ISR = niter,HSR = niter,h = niter)# for loop to go through each iteration and calculate the difference between # female and male survival rates for each stage.for(i in1:niter){ Adult <-pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="adult"& rate =="survival"& sex =="Female"), value) -pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="adult"& rate =="survival"& sex =="Male"), value) Fledgling <-pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="fledgling"& rate =="survival"& sex =="Female"), value) -pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="fledgling"& rate =="survival"& sex =="Male"), value) Groundling <-pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="groundling"& rate =="survival"& sex =="Female"), value) -pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="groundling"& rate =="survival"& sex =="Male"), value) Nestling <-pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="nestling"& rate =="survival"& sex =="Female"), value) -pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="nestling"& rate =="survival"& sex =="Male"), value) HSR <-0.5-pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="HSR"), value) h <-1-pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], stage =="h"), value) sex_diff_surv_output[i, 1] <-ifelse(length(Adult) >0, Adult, NA) sex_diff_surv_output[i, 2] <-ifelse(length(Fledgling) >0, Fledgling, NA) sex_diff_surv_output[i, 3] <-ifelse(length(Groundling) >0, Groundling, NA) sex_diff_surv_output[i, 4] <-ifelse(length(Nestling) >0, Nestling, NA) sex_diff_surv_output[i, 6] <-ifelse(length(HSR) >0, HSR, NA) sex_diff_surv_output[i, 7] <-ifelse(length(h) >0, h, NA) }# restructure the output and lable columns sex_diff_surv_output <- reshape2::melt(data = sex_diff_surv_output)colnames(sex_diff_surv_output) <-c("stage", "difference")# return the output sex_diff_surv_output}
Wrangling survival data
The following chunks show how we wrangled the raw field data to consolidate the data into the appropriate format for the survival analyses
Offspring encounter and census history
Code
status_dat_all <-# read raw dataread.csv("data/raw/Coucal_chick_survival_2001-2019_20200129.csv", header =TRUE, stringsAsFactors =FALSE, na.strings =c("", " ", "NA")) %>%# rename ring_ID column dplyr::rename(ring_ID = Ring_ID) %>%# make all entries lower case for consistencymutate(Fledged_status =tolower(Fledged.),site =tolower(site)) %>%# select variables of interest dplyr::select(species, ring_ID, lab_no, sex, year, site, nest_ID, pref_age, Fledged_status, postf_age, postf_status, ageC, statusC, lay_date, hatch_order) %>%# remove all white space from datamutate(across(everything(), ~str_trim(.x))) %>%mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%# specify empty data as NAmutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%# exclude all individuals that died in the nest# filter(Fledged_status == "yes") %>% # classify columnsmutate(sex =as.factor(sex),ageC =as.numeric(ageC),postf_age =as.numeric(postf_age),postf_status =as.numeric(postf_status),hatch_order =as.numeric(hatch_order),pref_age =as.numeric(pref_age)) %>%# remove rows with missing sex, age, and status infofilter(!is.na(sex) &!is.na(ageC) &!is.na(statusC)) %>%# make a unique id for each individualmutate(ind_ID =paste(nest_ID, lab_no, ring_ID, sep ="_"),# create the age of entry into the data (all at age 15)entry =0,# specify the age of death or censoringexit = ageC,# make the event numeric and specify if # the individual died (1) or was censored (0)event =as.numeric(statusC)) %>%# consolidate to variables of interest dplyr::select(species, ind_ID, nest_ID, year, sex, entry, exit, event, hatch_order)
Adult encounter history
Code
detect_dat <-read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na ="NA", col_types ="text") %>% dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>%mutate(month =str_sub(date_dec, start =5, end =6),day =str_sub(date_dec, start =7, end =8)) %>%mutate(date =as.Date(paste(year, month, day, sep ="-"), format ="%Y-%m-%d")) %>% dplyr::select(-month, -day, -date_dec) %>%mutate(across(c("sex", "site", "age_status"), tolower)) %>%mutate(sex =ifelse(sex =="female", "F", ifelse(sex =="male", "M", "XXX")),age_status =ifelse(age_status =="adult", "A", ifelse(age_status =="juvenile", "J", "XXX")),ring_ID =str_replace_all(string = ring_ID, fixed(" "), "")) %>%mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor))# make annual capture histories for adultsBC_detect_dat_A <- detect_dat %>%filter(age_status =="A"& species =="BC")# use the BaSTA function "CensusToCaptHist()" function to convert long format# encounter histories of each individual, to wide format with 1's and 0's for # each year of encounterBC_detect_dat_A_ch <-CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,d = BC_detect_dat_A$year) %>% dplyr::rename(ring_ID = ID) %>%mutate(ID =rownames(.)) %>%left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by ="ring_ID") %>%distinct()Black_Coucal_adult_CJS_ch <-data.frame(ch =apply(BC_detect_dat_A_ch[, 2:19 ] , 1, paste, collapse ="")) %>%bind_cols(., dplyr::select(BC_detect_dat_A_ch, sex, ring_ID)) %>%mutate(across(everything(), ~str_trim(.x))) %>%mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%distinct()WBC_detect_dat_A <- detect_dat %>%filter(age_status =="A"& species =="WBC") %>%mutate(year =as.numeric(year))WBC_detect_dat_A_ch <-CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,d = WBC_detect_dat_A$year) %>% dplyr::rename(ring_ID = ID) %>%mutate(ID =rownames(.)) %>%left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by ="ring_ID") %>%distinct()White_browed_Coucal_adult_CJS_ch <-data.frame(ch =apply(WBC_detect_dat_A_ch[, 2:16 ] , 1, paste, collapse ="")) %>%bind_cols(., dplyr::select(WBC_detect_dat_A_ch, sex, ring_ID)) %>%mutate(across(everything(), ~str_trim(.x))) %>%mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%distinct()
Bootstrapped mark-recapture analysis of White-browed Coucal adults
Here we resample the annual encounter history of adult White-browed Coucals and conduct a CJS mark-recapture analysis to derive the mean estimate and variance in annual apparent survival.
Set the ‘niter’ to run the bootstrap over several iterations (e.g., 10000). For reproducibility, set the seed with the set.seed() function. Set the ‘output.directory’ so that the subsequent chunks can refer to the appropriate folder within the project directory.
# save the outputsaveRDS(WBC_adult_survival_CJS_bootstrap_result,file ="output/bootstraps/raw/WBC_adult_survival_CJS_bootstrap_result.rds")# tidy it upWBC_adult_survival_CJS_bootstrap_tidy <-surv_boot_out_wrangle(species ="WBC", niter = niter, output_directory = output.directory)# save the tidy versionsaveRDS(WBC_adult_survival_CJS_bootstrap_tidy,file ="output/bootstraps/tidy/WBC_adult_survival_CJS_bootstrap_tidy.rds")
Derive the Black Coucal sex-specific survival rates by transforming White-browed Coucal sex-specific survival rates to reflect the ~0.7 ASR that is observed annually
Code
# bootstrapped WBC survival ratesWBC_adult_surival_boot_out <- WBC_boot_out$adult_reals_survival_rates_boot %>% dplyr::rename(value = estimate)# boostrapped BC survival rates (i.e., derived from the WBC rates)trans_WBC_adult_surival_boot_out <- WBC_boot_out$adult_reals_survival_rates_boot %>% dplyr::rename(value = estimate) %>%mutate(value =ifelse(sex =="Female", value *0.6, value *1.4))
Parameterizing demographic rates
Each of the following rates has a mean and standard deviation that are used in the stochastic simulation of adult sex ratio.
Mating system
Code
mating_dat <-# read raw dataread.csv("data/raw/Coucal_Nr_mates_2001_2019_20200123.csv", header =TRUE, stringsAsFactors =FALSE, na.strings =c("", " ", "NA")) %>%# make all entries lower case for consistencymutate(site =tolower(site),age_status =tolower(age_status),sex =tolower(sex)) %>%# specify empty data as NAmutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%# classify columnsmutate(sex =as.factor(sex),species =as.factor(species),Nr_partners =as.numeric(Nr_partners),age_status =as.factor(age_status),site =as.factor(site)) %>%# remove rows with missing sex and age, and post-fledged status infofilter(!is.na(sex) &!is.na(Nr_partners) &!is.na(species) & species !="CTC") %>%# make a unique id for each individualmutate(ind_ID =str_replace_all(ring_ID, fixed(" "), "")) %>%mutate(ID_partner.s. =str_replace_all(ID_partner.s., fixed(" "), "")) %>%# make a unique id for each individualmutate(ind_ID =str_replace_all(ind_ID, fixed("_"), "")) %>%mutate(ID_partner.s. =str_replace_all(ID_partner.s., fixed("_"), "")) %>%# consolidate to variables of interest dplyr::select(year, species, ind_ID, site, sex, age_status, Nr_partners, ID_partner.s.)sex_specific_mating_system <- mating_dat %>%group_by(ind_ID, sex, species) %>% dplyr::summarise(mean_annual_no_mates_ind =mean(Nr_partners, na.rm =TRUE),n_years =n_distinct(year)) %>%mutate(mu_i =ifelse(mean_annual_no_mates_ind <=1, 1, mean_annual_no_mates_ind)) %>%group_by(species, sex) %>% dplyr::summarise(mean_annual_no_mates =mean(mu_i, na.rm =TRUE),var_annual_no_mates =var(mu_i, na.rm =TRUE),median_annual_no_mates =median(mu_i, na.rm =TRUE),sd_annual_no_mates =sd(mu_i, na.rm =TRUE),n =n_distinct(ind_ID), .groups ="drop") %>%mutate(sample_size =paste("n = ", n, sep =""),species_plot =ifelse(species =="WBC", 1.9, 1.1),species_lab =ifelse(species =="WBC", 1, 2))# To obtain a female-based h index for each population the inverse of the mean # mu is calculated (i.e., Eq. 2 in manuscript)BC_h <-as.numeric(sex_specific_mating_system[which( sex_specific_mating_system$sex =="female"& sex_specific_mating_system$species =="BC"), "mean_annual_no_mates"])WBC_h <-as.numeric(sex_specific_mating_system[which( sex_specific_mating_system$sex =="female"& sex_specific_mating_system$species =="WBC"), "mean_annual_no_mates"])BC_h_sd <-as.numeric(sex_specific_mating_system[which( sex_specific_mating_system$sex =="female"& sex_specific_mating_system$species =="BC"), "sd_annual_no_mates"])WBC_h_sd <-as.numeric(sex_specific_mating_system[which( sex_specific_mating_system$sex =="female"& sex_specific_mating_system$species =="WBC"), "sd_annual_no_mates"])coucal_mating_system <-data.frame(trait =c("mating_system"),species =c("BC", "WBC"),mean =c(BC_h, WBC_h),sd =c(BC_h_sd, WBC_h_sd))mating_dat$species <-factor(mating_dat$species ,levels =c("BC","WBC"))mating_dat_plotting <- mating_dat %>%mutate(species_dummy =ifelse(species =="WBC", 1, 0)) %>%mutate(species_plot =ifelse(species =="WBC", 2.1, 0.9))mating_system_plot <- ggplot2::ggplot() +theme_bw() +geom_jitter(aes(y = Nr_partners, x = species_plot, color = sex),data = mating_dat_plotting, width =0.1, height =0.1,alpha =0.75, fill ="white", shape =16) +geom_pointrange(data = sex_specific_mating_system, aes(y = mean_annual_no_mates, x = species_plot, ymin = (mean_annual_no_mates-sd_annual_no_mates), ymax = (mean_annual_no_mates+sd_annual_no_mates),fill = sex),size =0.8, shape =21, color ="black") +geom_text(aes(y =0.2, x = species_lab, label = sample_size), data = sex_specific_mating_system,size =3, family ="Verdana") +geom_hline(yintercept =1.5, linetype ="dashed",alpha =0.5, color ="grey20") +facet_grid(. ~ sex, labeller =as_labeller(sex_names)) + luke_theme +theme(legend.position ="none",axis.title.x =element_blank(),strip.background =element_blank(),strip.text =element_text(size =12, face ="italic")) +scale_x_continuous(limits =c(0.5, 2.5), breaks =c(1, 2), labels = species_names) +scale_y_continuous(limits =c(0, 5.5), expand =c(0, 0), breaks =c(1, 2, 3, 4, 5),labels =c(1, 2, 3, 4, 5)) +ylab("Number of unique mates\nper year (± 1 SD)") +scale_color_manual(values = sex_pal2) +scale_fill_manual(values = sex_pal2) +annotate(geom ="text", y =1.25, x =1.5,label ="monogamy", family ="Verdana", color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =1.75, x =1.5,label ="polygamy", family ="Verdana",color ="black", size =3, fontface ='italic', hjust =0.5)
Clutch size
Code
# import raw csv data into Regg_data <-read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", na ="NA", col_types ="text") %>% dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,`unkn_sex_eggs+chiks`) %>%filter(species !="CTC") %>%filter(!is.na(nest_success)) %>% dplyr::rename(n_eggs =`Total clutch size`,n_chicks = all_chicks,male_eggs =`male_chi+eggs`, female_eggs =`female_chi+eggs`,unknown_eggs =`unkn_sex_eggs+chiks`) %>% dplyr::mutate(n_eggs =na_if(n_eggs, "?"),n_chicks =na_if(n_chicks, "?")) %>% dplyr::mutate(n_eggs =as.numeric(n_eggs),n_chicks =as.numeric(n_chicks),male_eggs =as.numeric(male_eggs),female_eggs =as.numeric(female_eggs),unknown_eggs =as.numeric(unknown_eggs)) %>%filter(n_eggs >= n_chicks)#### clutch size summary ----coucal_clutch_size <-data.frame(trait =c("clutch_size"),species =c("BC", "WBC"),mean =c(filter(egg_data, species =="BC") %>%summarise(mean_n_eggs =mean(n_eggs, na.rm =TRUE)) %>%pull(mean_n_eggs),filter(egg_data, species =="WBC") %>%summarise(mean_n_eggs =mean(n_eggs, na.rm =TRUE)) %>%pull(mean_n_eggs)),sd =c(filter(egg_data, species =="BC") %>%summarise(mean_n_eggs =sd(n_eggs, na.rm =TRUE)) %>%pull(mean_n_eggs),filter(egg_data, species =="WBC") %>%summarise(mean_n_eggs =sd(n_eggs, na.rm =TRUE)) %>%pull(mean_n_eggs)),n_nests =c(filter(egg_data, species =="BC") %>%summarise(n_nests =n_distinct(nest_ID)) %>%pull(n_nests),filter(egg_data, species =="WBC") %>%summarise(n_nests =n_distinct(nest_ID)) %>%pull(n_nests)),n_years =c(filter(egg_data, species =="BC") %>%summarise(n_nests =n_distinct(year)) %>%pull(n_nests),filter(egg_data, species =="WBC") %>%summarise(n_nests =n_distinct(year)) %>%pull(n_nests)))
Egg hatchability
Code
# import raw csv data into Regg_data <-read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", na ="NA", col_types ="text") %>% dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,`unkn_sex_eggs+chiks`) %>%filter(species !="CTC") %>%filter(!is.na(nest_success)) %>% dplyr::rename(n_eggs =`Total clutch size`,n_chicks = all_chicks,male_eggs =`male_chi+eggs`, female_eggs =`female_chi+eggs`,unknown_eggs =`unkn_sex_eggs+chiks`) %>% dplyr::mutate(n_eggs =na_if(n_eggs, "?"),n_chicks =na_if(n_chicks, "?")) %>% dplyr::mutate(n_eggs =as.numeric(n_eggs),n_chicks =as.numeric(n_chicks),male_eggs =as.numeric(male_eggs),female_eggs =as.numeric(female_eggs),unknown_eggs =as.numeric(unknown_eggs)) %>%filter(n_eggs >= n_chicks) %>%arrange(nest_ID)#### egg hatchability summary ----egg_surv_data <- egg_data %>% dplyr::mutate(n_dead_eggs = n_eggs - n_chicks,n_alive_eggs = n_chicks) %>%arrange(nest_ID)#### sample size summary ----# egg_surv_data %>% # dplyr::group_by(species) %>% # dplyr::summarise(n_nests = n_distinct(nest_ID))BC_egg_survival <- lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~1+ (1| nest_ID), data =filter(egg_surv_data, species =="BC"), family = binomial)WBC_egg_survival <- lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~1+ (1| nest_ID), data =filter(egg_surv_data, species =="WBC"), family = binomial)# # run tidy bootstrap to obtain model diagnosticstidy_BC_egg_survival <-tidy(BC_egg_survival, conf.int =TRUE, conf.method ="boot", nsim =1000)tidy_WBC_egg_survival <-tidy(WBC_egg_survival, conf.int =TRUE, conf.method ="boot", nsim =1000)coucal_egg_survival <-data.frame(trait =c("egg_survival"),species =c("BC", "WBC"),mean =c(invlogit(model_parameters(BC_egg_survival)$Coefficient)[1],invlogit(model_parameters(WBC_egg_survival)$Coefficient)[1]),CI_low =c(invlogit(model_parameters(BC_egg_survival)$CI_low)[1],invlogit(model_parameters(WBC_egg_survival)$CI_low)[1]),CI_high =c(invlogit(model_parameters(BC_egg_survival)$CI_high)[1],invlogit(model_parameters(WBC_egg_survival)$CI_high)[1]),n_nests =c(filter(egg_data, species =="BC") %>%summarise(n_nests =n_distinct(nest_ID)) %>%pull(n_nests),filter(egg_data, species =="WBC") %>%summarise(n_nests =n_distinct(nest_ID)) %>%pull(n_nests)),n_years =c(filter(egg_data, species =="BC") %>%summarise(n_nests =n_distinct(year)) %>%pull(n_nests),filter(egg_data, species =="WBC") %>%summarise(n_nests =n_distinct(year)) %>%pull(n_nests))) %>%mutate(sd =ifelse(!is.na(CI_low), approx_sd(x1 = CI_low, x2 = CI_high), CI_low))
#### Modeling (WBC) ----# "Fledge_age" as dependent variable, interaction with sex and speciesmod_fledge_age_WBC <-lmer(Fledge_age ~ sex + (1| nest_ID) + (1| year), data =filter(fledge_dat, species =="WBC"))# extract model coefficientsmod_fledge_age_WBC_coefs <-model_parameters(mod_fledge_age_WBC) %>%as.data.frame(.)# run tidy bootstrap to obtain model diagnosticstidy_mod_fledge_age_WBC <-tidy(mod_fledge_age_WBC, conf.int =TRUE, conf.method ="boot", nsim =1000)# run rptR to obtain repeatabilities of random effectsrpt_mod_fledge_age_WBC <-rpt(Fledge_age ~ sex + (1| nest_ID) + (1| year),grname =c("nest_ID", "year", "Fixed"),data =filter(fledge_dat, species =="WBC"),datatype ="Gaussian",nboot =1000, npermut =1000, ratio =TRUE,adjusted =TRUE, ncores =4, parallel =TRUE)# run partR2 to obtain marginal R2, parameter estimates, and beta weightsR2m_mod_fledge_age_WBC <-partR2(mod_fledge_age_WBC,partvars =c("sex"),R2_type ="marginal",nboot =1000, CI =0.95, max_level =1)R2c_mod_fledge_age_WBC <-partR2(mod_fledge_age_WBC,partvars =c("sex"),R2_type ="conditional",nboot =1000, CI =0.95, max_level =1)# save model, tidy, rptR, and partR2 output as a liststats_mod_fledge_age_WBC <-list(mod = mod_fledge_age_WBC,tidy = tidy_mod_fledge_age_WBC,rptR = rpt_mod_fledge_age_WBC,partR2m = R2m_mod_fledge_age_WBC,partR2c = R2c_mod_fledge_age_WBC)
pull species- and sex-specific means and sds for each parameter
Code
k_dist_BC =c(pull(filter(coucal_clutch_size, species =="BC"), mean), pull(filter(coucal_clutch_size, species =="BC"), sd))HSR_dist_BC =c(pull(filter(coucal_HSR, species =="BC"), mean), pull(filter(coucal_HSR, species =="BC"), sd))h_dist_BC =c(pull(filter(coucal_mating_system, species =="BC"), mean), pull(filter(coucal_mating_system, species =="BC"), sd))egg_surv_dist_BC =c(pull(filter(coucal_egg_survival, species =="BC"), mean),pull(filter(coucal_egg_survival, species =="BC"), sd))fledge_age_distF_BC =c(pull(filter(coucal_fledge_age, species =="BC"& sex =="F"), mean), pull(filter(coucal_fledge_age, species =="BC"& sex =="F"), sd))fledge_age_distM_BC =c(pull(filter(coucal_fledge_age, species =="BC"& sex =="M"), mean), pull(filter(coucal_fledge_age, species =="BC"& sex =="M"), sd))flight_age_distF_BC =c(pull(filter(coucal_flight_age, species =="BC"& sex =="F"), mean), pull(filter(coucal_flight_age, species =="BC"& sex =="F"), sd))flight_age_distM_BC =c(pull(filter(coucal_flight_age, species =="BC"& sex =="M"), mean), pull(filter(coucal_flight_age, species =="BC"& sex =="M"), sd))
pull species- and sex-specific means and sds for each parameter
Code
k_dist_WBC =c(pull(filter(coucal_clutch_size, species =="WBC"), mean), pull(filter(coucal_clutch_size, species =="WBC"), sd))HSR_dist_WBC =c(pull(filter(coucal_HSR, species =="WBC"), mean), pull(filter(coucal_HSR, species =="WBC"), sd))h_dist_WBC =c(pull(filter(coucal_mating_system, species =="WBC"), mean), pull(filter(coucal_mating_system, species =="WBC"), sd))egg_surv_dist_WBC =c(pull(filter(coucal_egg_survival, species =="WBC"), mean),pull(filter(coucal_egg_survival, species =="WBC"), sd))fledge_age_distF_WBC =c(pull(filter(coucal_fledge_age, species =="WBC"& sex =="F"), mean), pull(filter(coucal_fledge_age, species =="WBC"& sex =="F"), sd))fledge_age_distM_WBC =c(pull(filter(coucal_fledge_age, species =="WBC"& sex =="M"), mean), pull(filter(coucal_fledge_age, species =="WBC"& sex =="M"), sd))flight_age_distF_WBC =c(pull(filter(coucal_flight_age, species =="WBC"& sex =="F"), mean), pull(filter(coucal_flight_age, species =="WBC"& sex =="F"), sd))flight_age_distM_WBC =c(pull(filter(coucal_flight_age, species =="WBC"& sex =="M"), mean), pull(filter(coucal_flight_age, species =="WBC"& sex =="M"), sd))
Visualize the sex-specific variation in cumulative survival and hazard rate of premature individuals derived from the stochastic simulation
Code
grp_means <-bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>%filter(rate =="development") %>%pivot_wider(names_from = stage, values_from = value) %>%mutate(flight_age = fledge_age + flight_age) %>%pivot_longer(c(fledge_age, flight_age), names_to ="stage", values_to ="value") %>%mutate(stage_sex =paste(sex, stage, sep ="_")) %>%group_by(species, stage_sex, sex, stage) %>% dplyr::summarise(grp.mean =mean(value))BC_density_plotA <-bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>%filter(rate =="development"& species =="BC") %>%pivot_wider(names_from = stage, values_from = value) %>%mutate(flight_age = fledge_age + flight_age) %>%pivot_longer(c(fledge_age, flight_age), names_to ="stage", values_to ="value") %>%mutate(stage_sex =paste(sex, stage, sep ="_")) %>%ggplot() +geom_density(aes(value, y=..scaled.., fill =fct_rev(stage_sex)),alpha =0.7, color ="grey40", size =0.3) +theme_void() +theme(legend.position ="none",plot.margin =unit(c(0,0,0,0.5), "cm"),plot.background =element_rect(fill ='transparent', color =NA)) +scale_fill_manual(values =rev(sex_pal3)) +scale_x_continuous(limits =c(0, 70),expand =c(0.01, 0.01)) +ggtitle('A')BC_density_plotB <-bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>%filter(rate =="development"& species =="BC") %>%pivot_wider(names_from = stage, values_from = value) %>%mutate(flight_age = fledge_age + flight_age) %>%pivot_longer(c(fledge_age, flight_age), names_to ="stage", values_to ="value") %>%mutate(stage_sex =paste(sex, stage, sep ="_")) %>%ggplot() +geom_density(aes(value, y=..scaled.., fill =fct_rev(stage_sex)),alpha =0.7, color ="grey40", size =0.3) +theme_void() +theme(legend.position ="none",plot.margin =unit(c(0,0,0,0.5), "cm"),plot.background =element_rect(fill ='transparent', color =NA)) +scale_fill_manual(values =rev(sex_pal3)) +scale_x_continuous(limits =c(0, 70),expand =c(0.01, 0.01)) +ggtitle('B')WBC_density_plot <-bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>%filter(rate =="development"& species =="WBC") %>%pivot_wider(names_from = stage, values_from = value) %>%mutate(flight_age = fledge_age + flight_age) %>%pivot_longer(c(fledge_age, flight_age), names_to ="stage", values_to ="value") %>%mutate(stage_sex =paste(sex, stage, sep ="_")) %>%ggplot() +geom_density(aes(value, y=..scaled.., fill =fct_rev(stage_sex)),alpha =0.7, color ="grey40", size =0.3) +theme_void() +theme(legend.position ="none",plot.margin =unit(c(0,0,0,0.5), "cm"),plot.background =element_rect(fill ='transparent', color =NA)) +scale_fill_manual(values =rev(sex_pal3)) +scale_x_continuous(limits =c(0, 70),expand =c(0.01, 0.03))age_means_BC <- BC_hazard_rate_boot_tidy$hazard_rates_boot %>%filter(iter %in%c(as.character(seq(from =1, to =1000, by =1)))) %>% dplyr::group_by(age, sex) %>% dplyr::summarise(mean_rate =mean(estimate, na.rm =TRUE))age_means_WBC <- WBC_hazard_rate_boot_tidy$hazard_rates_boot %>%filter(iter %in%c(as.character(seq(from =1, to =1000, by =1)))) %>% dplyr::group_by(age, sex) %>% dplyr::summarise(mean_rate =mean(estimate, na.rm =TRUE))BC_surv_plot <-ggplot() + luke_theme +theme(legend.position ="none",strip.background =element_blank(),strip.text =element_text(size =12, face ="italic"),axis.text.x =element_blank(),axis.title.x =element_blank(),axis.title.y =element_blank(),plot.margin =unit(c(0,0,0.5,0.5), "cm"),axis.ticks =element_line(size =0.25, colour ="grey40"),axis.ticks.length =unit(0.1, "cm"),axis.text.y =element_text(size =8),panel.grid.major.x =element_line(colour ="grey70", size =0.1),plot.background =element_rect(fill ='transparent', color =NA)) +geom_vline(data =filter(grp_means, species =="BC"& stage =="fledge_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1) +geom_vline(data =filter(grp_means, species =="BC"& stage =="flight_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1) +geom_line(data =filter(BC_hazard_rate_boot_tidy$hazard_rates_boot, iter %in%c(as.character(seq(from =1, to =1000, by =1)))),aes(x = age, y = estimate, group =interaction(iter, sex), color = sex),alpha =0.05, size =0.25) +geom_line(data = age_means_BC,aes(x = age, y = mean_rate, group = sex), color =c(rep("#00496C", 70), rep("#E5BD3A", 70)),alpha =1) +facet_grid(species ~ ., labeller =as_labeller(species_names)) +scale_colour_manual(values = sex_pal2) +ylab("Estimated daily survival rate") +xlab("Age (Days since hatching)") +scale_y_continuous(limits =c(0.9, 1),breaks =seq(from =0.9, to =1, by =0.025)) +scale_x_continuous(limits =c(0, 70), breaks =seq(0, 70, by =5), expand =c(0.01, 0.01),labels =as.character(seq(0, 70, by =5))) +annotate(geom ="text", y =0.91, x =mean(c(pull(filter(grp_means, species =="BC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="BC"& stage =="fledge_age")[2,5], grp.mean)))/2,label ="nestling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.91, x =mean(c(pull(filter(grp_means, species =="BC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="BC"& stage =="fledge_age")[2,5], grp.mean))) + ((mean(c(pull(filter(grp_means, species =="BC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="BC"& stage =="flight_age")[2,5], grp.mean))) -mean(c(pull(filter(grp_means, species =="BC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="BC"& stage =="fledge_age")[2,5], grp.mean))))/2),label ="groundling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.91, x =mean(c(pull(filter(grp_means, species =="BC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="BC"& stage =="flight_age")[2,5], grp.mean))) + ((70-mean(c(pull(filter(grp_means, species =="BC"& stage =="flight_age")[1,5], grp.mean), pull(filter(grp_means, species =="BC"& stage =="flight_age")[2,5], grp.mean))))/2),label ="fledgling",color ="black", size =3, fontface ='italic', hjust =0.5)WBC_surv_plot <-ggplot() + luke_theme +theme(legend.position ="none",strip.background =element_blank(),strip.text =element_text(size =12, face ="italic"),plot.margin =unit(c(0,0,0.5,0.5), "cm"),axis.ticks =element_line(size =0.25, colour ="grey40"),axis.ticks.length =unit(0.1, "cm"),axis.text.x =element_text(size =8),axis.text.y =element_text(size =8),panel.grid.major.x =element_line(colour ="grey70", size =0.1),plot.background =element_rect(fill ='transparent', color =NA)) +geom_vline(data =filter(grp_means, species =="WBC"& stage =="fledge_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1) +geom_vline(data =filter(grp_means, species =="WBC"& stage =="flight_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1) +geom_line(data =filter(WBC_hazard_rate_boot_tidy$hazard_rates_boot, iter %in%c(as.character(seq(from =1, to =1000, by =1)))),aes(x = age, y = estimate, group =interaction(iter, sex), color = sex),alpha =0.05, size =0.25) +geom_line(data = age_means_WBC,aes(x = age, y = mean_rate, group = sex), color =c(rep("#00496C", 70), rep("#E5BD3A", 70)),alpha =1) +facet_grid(species ~ ., labeller =as_labeller(species_names)) +scale_colour_manual(values = sex_pal2) +ylab(" Estimated daily survival rate") +xlab("Age (days since hatching)") +scale_y_continuous(limits =c(0.9, 1),breaks =seq(from =0.9, to =1, by =0.025)) +scale_x_continuous(limits =c(0, 70), breaks =seq(0, 70, by =5), expand =c(0.01, 0.01),labels =as.character(seq(0, 70, by =5))) +annotate(geom ="text", y =0.91, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean)))/2,label ="nestling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.91, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean))) + ((mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))) -mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean))))/2),label ="groundling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.91, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))) + ((70-mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean), pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))))/2),label ="fledgling",color ="black", size =3, fontface ='italic', hjust =0.5)# BC_coxme <- coxme(Surv(exit, event) ~ sex + (1|year) + (1|nest_ID), data = status_dat_all %>% filter(species == "BC")) # summary(BC_coxme)BC_survfit_sex <-survfit(Surv(exit, event) ~ sex , data = status_dat_all %>%filter(species =="BC"))BC_plot <-ggsurvplot(BC_survfit_sex,risk.table = F,ncensor.plot = F,conf.int = T,ylab ="Cumulative Survival",xlab ="Age (Days since hatching)",xlim =c(0,70),break.time.by =10,palette = sex_pal2,ggtheme = luke_theme +theme(plot.title =element_text(hjust =0.5, face ="italic"),legend.background =element_rect(fill="transparent")),legend =c(0.7, 0.9),font.legend =c(10, "plain", "black"),legend.title ="",legend.labs =c("Female", "Male"))# WBC_coxme <- coxme(Surv(exit, event) ~ sex + (1|year) + (1|nest_ID), data = status_dat_all %>% filter(species == "WBC")) # summary(WBC_coxme)WBC_survfit_sex <-survfit(Surv(exit, event) ~ sex , data = status_dat_all %>%filter(species =="WBC"))WBC_plot <-ggsurvplot(WBC_survfit_sex,risk.table = F,ncensor.plot = F,conf.int = T,ylab ="Cumulative Survival",xlab ="Age (Days since hatching)",xlim =c(0,70),break.time.by =10,palette = sex_pal2,ggtheme = luke_theme +theme(plot.title =element_text(hjust =0.5, face ="italic"),legend.background =element_rect(fill="transparent")),legend =c(0.7, 0.9),font.legend =c(10, "plain", "black"),legend.title ="",legend.labs =c("Female", "Male"))BC_plot1 <- BC_plot$plot +geom_vline(data =filter(grp_means, species =="BC"& stage =="fledge_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1, color =c("#be9c2e", "#016392")) +geom_vline(data =filter(grp_means, species =="BC"& stage =="flight_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1, color =c("#be9c2e", "#016392")) + luke_theme +theme(legend.position =c(0.8, 0.8),plot.margin =unit(c(0,0,0.5,0), "cm"),legend.text=element_text(size =10),legend.title =element_blank(),legend.key.size =unit(0.3, 'cm'),# panel.border = element_blank(),strip.background =element_blank(),strip.text =element_text(size =12, face ="italic"),axis.title.x =element_blank(),axis.title.y =element_blank(),# axis.ticks.y = element_blank(),axis.ticks =element_line(size =0.25, colour ="grey40"),axis.ticks.length =unit(0.1, "cm"),axis.text.x =element_blank(),axis.text.y =element_text(size =8),panel.grid.major.x =element_line(colour ="grey70", size =0.1),# plot.title = element_text(hjust = 0.5, face = "italic"),legend.background =element_rect(fill="transparent"),plot.background =element_rect(fill ='transparent', color =NA)) +scale_x_continuous(limits =c(0, 70), breaks =seq(0, 70, by =5), expand =c(0.01, 0.01),labels =as.character(seq(0, 70, by =5))) +#c(-60, -30, 0, 30, 60)) +# ylab(expression(paste("Cumulative survival" %+-% "95% CI", sep = ""))) +# xlab("Age (Days since hatching)") +scale_color_manual(values = sex_pal2,labels =c("Female", "Male")) +scale_fill_manual(values = sex_pal2,labels =c("Female", "Male")) +annotate(geom ="text", y =0.08, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean)))/2,label ="nestling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.08, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean))) + ((mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))) -mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean))))/2),label ="groundling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.08, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))) + ((70-mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean), pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))))/2),label ="fledgling",color ="black", size =3, fontface ='italic', hjust =0.5)WBC_plot1 <- WBC_plot$plot +geom_vline(data =filter(grp_means, species =="WBC"& stage =="fledge_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1, color =c("#be9c2e", "#016392")) +geom_vline(data =filter(grp_means, species =="WBC"& stage =="flight_age"),aes(xintercept = grp.mean, color = sex), linetype ="dashed", alpha =1, color =c("#be9c2e", "#016392")) + luke_theme +theme(legend.position ="none",plot.margin =unit(c(0,0,0.5,0), "cm"),# panel.border = element_blank(),strip.background =element_blank(),strip.text =element_text(size =12, face ="italic"),axis.title.x =element_text(size =12),axis.ticks =element_line(size =0.25, colour ="grey40"),axis.ticks.length =unit(0.1, "cm"),axis.text.x =element_text(size =8),axis.text.y =element_text(size =8),# axis.ticks.y = element_blank(),legend.background =element_blank(),panel.grid.major.x =element_line(colour ="grey70", size =0.1),plot.background =element_rect(fill ='transparent', color =NA)) +# axis.title.y = element_blank(),# axis.text.y = element_blank()) +scale_x_continuous(limits =c(0, 70), breaks =seq(0, 70, by =5), expand =c(0.01, 0.01),labels =as.character(seq(0, 70, by =5))) +#c(-60, -30, 0, 30, 60)) +# xlab("Age (Days since hatching)") +scale_color_manual(values = sex_pal2) +scale_fill_manual(values = sex_pal2) +ylab(" Cumulative survival ± 95% CI") +xlab("Age (days since hatching)") +annotate(geom ="text", y =0.08, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean)))/2,label ="nestling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.08, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean))) + ((mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))) -mean(c(pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="fledge_age")[2,5], grp.mean))))/2),label ="groundling",color ="black", size =3, fontface ='italic', hjust =0.5) +annotate(geom ="text", y =0.08, x =mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean),pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))) + ((70-mean(c(pull(filter(grp_means, species =="WBC"& stage =="flight_age")[1,5], grp.mean), pull(filter(grp_means, species =="WBC"& stage =="flight_age")[2,5], grp.mean))))/2),label ="fledgling",color ="black", size =3, fontface ='italic', hjust =0.5)
Cox-proportional hazards model
Here we statistically test for sex-specific age-dependent mortality hazard, while accounting for non-independence among siblings by clustering by nest_ID. We also test the assumptions of proportional hazards based on scaled Schoenfeld residuals (met: black coucal: P = 0.5; white-browed coucal: P = 0.08). In summary, we find no sex effect on hazard rate in either species (black coucal: eβ = 0.91, 95% CI = 0.70–1.19, P = 0.49; white-browed coucal: eβ = 0.90, 95% CI = 0.67–1.19, P = 0.47).
Code
BC_cox_ph_sex <-coxph(Surv(time = entry, time2 = exit, event = event) ~ sex +cluster(nest_ID), data = status_dat_all %>%filter(species =="BC"))WBC_cox_ph_sex <-coxph(Surv(time = entry, time2 = exit, event = event) ~ sex +cluster(nest_ID), data = status_dat_all %>%filter(species =="WBC"))
Black Coucal offspring
no sex effect on hazard rate: eβ = 0.91, 95% CI = 0.70–1.19, P = 0.493
Code
summary(BC_cox_ph_sex)
Call:
coxph(formula = Surv(time = entry, time2 = exit, event = event) ~
sex, data = status_dat_all %>% filter(species == "BC"), cluster = nest_ID)
n= 636, number of events= 208
coef exp(coef) se(coef) robust se z Pr(>|z|)
sexM -0.09308 0.91112 0.13911 0.13589 -0.685 0.493
exp(coef) exp(-coef) lower .95 upper .95
sexM 0.9111 1.098 0.6981 1.189
Concordance= 0.503 (se = 0.019 )
Likelihood ratio test= 0.45 on 1 df, p=0.5
Wald test = 0.47 on 1 df, p=0.5
Score (logrank) test = 0.45 on 1 df, p=0.5, Robust = 0.47 p=0.5
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).
Schoenfeld residuals and diagnostic plot show that assumptions are met (P = 0.5)
Survival dynamics of radio-marked Black (top) and White-browed (bottom) Coucal offspring over the first 70 days since hatching. A) Kaplan-Meier plots of the sex-specific cumulative survival. B) Inverted sex-specific smoothed hazard functions (right) showing the estimated daily survival rate. Fine lines show the estimated function of survival for a bootstrap that randomly sampled one offspring per brood to account for non-independence. Thick lines show the average trend for the bootstrapping procedure. Density plots show the sex-specific posterior distributions of age-points of nest departure and first flight given the bootstrap simulation.
Visualize the sex- and stage-specific variation in the survival rates derived from the stochastic simulation
Sensitivity analysis function takes the vital rate summary of the bootstrap procedure and conducts a perturbation analysis on each rate to assess how a proportional change in a given vital rate changes the ASR
arguments:
vital_rate_summary: object (a list of all vital rates) produced by the make_mprime_matrix() and make_treat_matrix() functions below
matrix_str: the matrix_structure expression
h: the mating system index specified within vital_rate_summary
k: the clutch size specified within vital_rate_summary
HSR: the hatching sex ratio specified within vital_rate_summary
niter: number of time steps to project matrix model
ASR: the stable stage distribution ASR derived from the matrix_ASR() function applied to the mprime or treat matrix
lambda: the stable stage distribution lambda derived from the matrix_ASR() function applied to the mprime or treat matrix
Code
sensitivity_analysis <-function(vital_rate_summary, matrix_str, h =1, k =4, HSR, niter =1000, ASR, lambda){# make a list of all parameters vr <-list(F_Nestling_survival = vital_rate_summary$F_Nestling_survival,F_Groundling_survival = vital_rate_summary$F_Groundling_survival,F_Fledgling_survival = vital_rate_summary$F_Fledgling_survival,F_Adult_survival = vital_rate_summary$F_Adult_survival,M_Nestling_survival = vital_rate_summary$M_Nestling_survival,M_Groundling_survival = vital_rate_summary$M_Groundling_survival,M_Fledgling_survival = vital_rate_summary$M_Fledgling_survival,M_Adult_survival = vital_rate_summary$M_Adult_survival)# number of stages in the matrix no_stages <-sqrt(length(matrix_str))# Define plover life-stages of the Ceuta snowy plover matrix model stages <-c("F_1st_year", "F_Adult", "M_1st_year", "M_Adult")# an empty t by x matrix stage <-matrix(numeric(no_stages * niter), nrow = no_stages)# an empty t vector to store the population sizes pop <-numeric(niter)# dataframe to store the perturbation results ASR_pert_results <-data.frame(parameter =c("F_Nestling_survival", "F_Groundling_survival", "F_Fledgling_survival", "F_Adult_survival","M_Nestling_survival", "F_Groundling_survival", "M_Fledgling_survival", "M_Adult_survival","h", "k", "HSR", "ISR"),sensitivities =numeric(length(vr) +4),elasticities =numeric(length(vr) +4)) lambda_pert_results <-data.frame(parameter =c("F_Nestling_survival", "F_Groundling_survival", "F_Fledgling_survival", "F_Adult_survival","M_Nestling_survival", "F_Groundling_survival", "M_Fledgling_survival", "M_Adult_survival","h", "k", "HSR", "ISR"),sensitivities =numeric(length(vr) +4),elasticities =numeric(length(vr) +4))# specifiy how many survival rates there are n <-length(vr)# create vectors of perturbations to test on parameters of the matrix model vr_nums <-seq(0, 1, 0.01) # proportional changes in survival and HSR (i.e., between 0 an 1) h_nums <-seq(0, 2, 0.02) # proportional changes in h index (i.e., between 0 and 2) k_nums <-seq(3, 5, 0.02) # proportional changes in k (i.e, between 3 and 5)# create empty dataframes to store the perturbation results for ASR and lambda vr_pert_ASR <-matrix(numeric(n *length(vr_nums)),ncol = n, dimnames =list(vr_nums, names(vr))) h_pert_ASR <-matrix(numeric(length(h_nums)),ncol =1, dimnames =list(h_nums, "h")) k_pert_ASR <-matrix(numeric(length(k_nums)),ncol =1, dimnames =list(k_nums, "k")) HSR_pert_ASR <-matrix(numeric(length(vr_nums)),ncol =1, dimnames =list(vr_nums, "HSR")) vr_pert_lambda <-matrix(numeric(n *length(vr_nums)),ncol = n, dimnames =list(vr_nums, names(vr))) h_pert_lambda <-matrix(numeric(length(h_nums)),ncol =1, dimnames =list(h_nums, "h")) k_pert_lambda <-matrix(numeric(length(k_nums)),ncol =1, dimnames =list(k_nums, "k")) HSR_pert_lambda <-matrix(numeric(length(vr_nums)),ncol =1, dimnames =list(vr_nums, "HSR"))#### perturbation of survival rates ####for (g in1:n) { # pick a column (i.e., a variable) vr2 <- vr # reset the vital rates to the originalfor (i in1:length(vr_nums)) # pick a perturbation level { vr2[[g]] <- vr_nums[i] # specify the vital rate with the new perturbation level A <-matrix(sapply(matrix_str, eval, vr2, NULL), nrow =sqrt(length(matrix_str)), byrow=TRUE, dimnames =list(stages, stages)) # build the matrix with the new value# reset the starting stage distribution for simulation (all with 10 individuals) m <-rep(10, no_stages) for (j in1:niter) { # project the matrix through t iteration# stage distribution at time t stage[,j] <- m# population size at time t pop[j] <-sum(m)# number of male adults at time t M2 <- stage[4, j] * A["M_Adult", "M_Adult"]# number of female adults at time t F2 <- stage[2, j] * A["F_Adult", "F_Adult"]# Female freq-dep fecundity of Female chicks A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1- HSR) )# Female freq-dep fecundity of Male chicks A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)# Male freq-dep fecundity of Female chicks A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1- HSR))# Male freq-dep fecundity of Male chicks A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)# define the new n (i.e., new stage distribution at time t) m <- A %*% (m) }# define rownames of stage matrixrownames(stage) <-rownames(A)# define colnames of stage matrixcolnames(stage) <-0:(niter -1)# calculate the proportional stable stage distribution stage <-apply(stage, 2, function(x) x /sum(x))# define stable stage as the last stage stable.stage <- stage[, niter]# calc ASR as the proportion of the adult stable stage class that is male vr_pert_ASR[i, g] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + stable.stage[no_stages])# calc lambda as the pop change in the counts of the last two iterations vr_pert_lambda[i, g] <- pop[niter]/pop[niter -1] }# get the spline function of ASR spl_ASR <-smooth.spline(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g] ~names(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g]))# estimate the slope of the tangent of the spline at the vital rate ASR_pert_results[g, 2] <-predict(spl_ASR, x = vr[[g]], deriv =1)$y# re-scale sensitivity into elasticity ASR_pert_results[g, 3] <- vr[[g]] / ASR * ASR_pert_results[g, 2]# do the same steps but for lambda spl_lambda <-smooth.spline(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g] ~names(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g])) lambda_pert_results[g, 2] <-predict(spl_lambda, x = vr[[g]], deriv =1)$y lambda_pert_results[g, 3] <- vr[[g]] / lambda * lambda_pert_results[g, 2] }#### perturbation of the h index parameter ####for (i in1:length(h_nums)) # pick a perturbation level { A <-matrix(sapply(matrix_str, eval, vr, NULL), nrow =sqrt(length(matrix_str)), byrow=TRUE, dimnames =list(stages, stages)) # reset the starting stage distribution for simulation (all with 10 individuals) m <-rep(10, no_stages) for (j in1:niter) { # project the matrix through t iteration# stage distribution at time t stage[,j] <- m# population size at time t pop[j] <-sum(m)# number of male adults at time t M2 <- stage[4, j] * A["M_Adult", "M_Adult"]# number of female adults at time t F2 <- stage[2, j] * A["F_Adult", "F_Adult"]# Female freq-dep fecundity of Female chicks A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * (1- HSR) )# Female freq-dep fecundity of Male chicks A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * HSR)# Male freq-dep fecundity of Female chicks A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * (1- HSR))# Male freq-dep fecundity of Male chicks A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * HSR)# define the new n (i.e., new stage distribution at time t) m <- A %*% (m) }# define rownames of stage matrixrownames(stage) <-rownames(A)# define colnames of stage matrixcolnames(stage) <-0:(niter -1)# calculate the proportional stable stage distribution stage <-apply(stage, 2, function(x) x /sum(x))# define stable stage as the last stage stable.stage <- stage[, niter]# calc ASR as the proportion of the adult stable stage class that is male h_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages /2] + stable.stage[no_stages])# calc lambda as the pop change in the counts of the last two iterations h_pert_lambda[i, ] <- pop[niter]/pop[niter -1] }# get the spline function of ASR spl_ASR <-smooth.spline(h_pert_ASR[which(!is.na(h_pert_ASR)), 1] ~names(h_pert_ASR[which(!is.na(h_pert_ASR)), ]))# estimate the slope of the tangent of the spline at the vital rate ASR_pert_results[n +1, 2] <-predict(spl_ASR, x = h, deriv =1)$y# re-scale sensitivity into elasticity ASR_pert_results[n +1, 3] <- h / ASR * ASR_pert_results[n +1, 2]# do the same steps but for lambda spl_lambda <-smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~names(h_pert_lambda[which(!is.na(h_pert_lambda)), ])) lambda_pert_results[n +1, 2] <-predict(spl_lambda, x = h, deriv =1)$y lambda_pert_results[n +1, 3] <- h / lambda * lambda_pert_results[n +1, 2]#### perturbation of k parameter ####for (i in1:length(k_nums)) # pick a perturbation level { A <-matrix(sapply(matrix_str, eval, vr, NULL), nrow =sqrt(length(matrix_str)), byrow=TRUE, dimnames =list(stages, stages))# reset the starting stage distribution for simulation (all with 10 individuals) m <-rep(10, no_stages) for (j in1:niter) { # project the matrix through t iteration# stage distribution at time t stage[,j] <- m# population size at time t pop[j] <-sum(m)# number of male adults at time t M2 <- stage[4, j] * A["M_Adult", "M_Adult"]# number of female adults at time t F2 <- stage[2, j] * A["F_Adult", "F_Adult"]# Female freq-dep fecundity of Female chicks A["F_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * (1- HSR) )# Female freq-dep fecundity of Male chicks A["M_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * HSR)# Male freq-dep fecundity of Female chicks A["F_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * (1- HSR))# Male freq-dep fecundity of Male chicks A["M_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * HSR)# define the new n (i.e., new stage distribution at time t) m <- A %*% (m) }# define rownames of stage matrixrownames(stage) <-rownames(A)# define colnames of stage matrixcolnames(stage) <-0:(niter -1)# calculate the proportional stable stage distribution stage <-apply(stage, 2, function(x) x /sum(x))# define stable stage as the last stage stable.stage <- stage[, niter]# calc ASR as the proportion of the adult stable stage class that is male k_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + stable.stage[no_stages])# calc lambda as the pop change in the counts of the last two iterations k_pert_lambda[i, ] <- pop[niter]/pop[niter -1] }# get the spline function of ASR spl_ASR <-smooth.spline(k_pert_ASR[which(!is.na(k_pert_ASR)), 1] ~names(k_pert_ASR[which(!is.na(k_pert_ASR)), ]))# estimate the slope of the tangent of the spline at the vital rate ASR_pert_results[n +2, 2] <-predict(spl_ASR, x = k, deriv =1)$y# re-scale sensitivity into elasticity ASR_pert_results[n +2, 3] <- k / ASR * ASR_pert_results[n +2, 2]# do the same steps but for lambda spl_lambda <-smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~names(h_pert_lambda[which(!is.na(h_pert_lambda)), ])) lambda_pert_results[n +2, 2] <-predict(spl_lambda, x = k, deriv =1)$y lambda_pert_results[n +2, 3] <- k / lambda * lambda_pert_results[n +2, 2]#### perturbation of HSR ####for (i in1:length(vr_nums)) # pick a perturbation level { A <-matrix(sapply(matrix_str, eval, vr, NULL), nrow =sqrt(length(matrix_str)), byrow=TRUE, dimnames =list(stages, stages))# reset the starting stage distribution for simulation (all with 10 individuals) m <-rep(10, no_stages) for (j in1:niter) { # project the matrix through t iteration# stage distribution at time t stage[,j] <- m# population size at time t pop[j] <-sum(m)# number of male adults at time t M2 <- stage[4, j] * A["M_Adult", "M_Adult"]# number of female adults at time t F2 <- stage[2, j] * A["F_Adult", "F_Adult"]# Female freq-dep fecundity of Female chicks A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1- vr_nums[i]) )# Female freq-dep fecundity of Male chicks A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * vr_nums[i])# Male freq-dep fecundity of Female chicks A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1- vr_nums[i]))# Male freq-dep fecundity of Male chicks A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * vr_nums[i])# define the new n (i.e., new stage distribution at time t) m <- A %*% (m) }# define rownames of stage matrixrownames(stage) <-rownames(A)# define colnames of stage matrixcolnames(stage) <-0:(niter -1)# calculate the proportional stable stage distribution stage <-apply(stage, 2, function(x) x /sum(x))# define stable stage as the last stage stable.stage <- stage[, niter]# calc ASR as the proportion of the adult stable stage class that is male HSR_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + stable.stage[no_stages])# calc lambda as the pop change in the counts of the last two iterations HSR_pert_lambda[i, ] <- pop[niter] / pop[niter -1] }# get the spline function of ASR spl_ASR <-smooth.spline(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), 1] ~names(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), ]))# estimate the slope of the tangent of the spline at the vital rate ASR_pert_results[n +3, 2] <-predict(spl_ASR, x = HSR, deriv =1)$y# re-scale sensitivity into elasticity ASR_pert_results[n +3, 3] <- HSR / ASR * ASR_pert_results[n +3, 2]# do the same steps but for lambda spl_lambda <-smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~names(h_pert_lambda[which(!is.na(h_pert_lambda)), ])) lambda_pert_results[n +3, 2] <-predict(spl_lambda, x = HSR, deriv =1)$y lambda_pert_results[n +3, 3] <- HSR/lambda * lambda_pert_results[n +3, 2]#### store all results into a list ---- result <-list(ASR_pert_results = ASR_pert_results,lambda_pert_results = lambda_pert_results,HSR_pert_ASR = HSR_pert_ASR,HSR_pert_lambda = HSR_pert_lambda,k_pert_ASR = k_pert_ASR,k_pert_lambda = k_pert_lambda,h_pert_ASR = h_pert_ASR,h_pert_lambda = h_pert_lambda,vr_pert_ASR = vr_pert_ASR,vr_pert_lambda = vr_pert_lambda) }
matrix_structure
matrix_structure is an expression that determines the structure of the two-sex matrix from a set of vital rates
Code
matrix_structure <-expression(# top row of matrix0, NA, 0, NA,# second row of matrix (F_Nestling_survival * F_Groundling_survival * F_Fledgling_survival), F_Adult_survival, 0, 0,# third row of matrix0, NA, 0, NA,# fourth row of matrix0, 0, (M_Nestling_survival * M_Groundling_survival * M_Fledgling_survival), M_Adult_survival)
make_treat_matrix()
This function makes a treatment matrix from a summary of the bootstrapped vital rates
arguments:
survival_rates_boot_summary: the $vital_rate_ests_boot element in the tidied _hazard_rate_boot_tidy object from above (a summary of the bootstrapped vital rates
species: species name (either “BC” or “WBC”)
h: the mating system index specified in the parameter_distributions
k: the clutch size specified in the parameter_distributions
HSR: the hatching sex ratio specified in the parameter_distributions
Code
make_treat_matrix <-function(survival_rates_boot_summary, species_name, h, k, HSR){list(F_Nestling_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_nestling_survival")[, 4],F_Groundling_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_groundling_survival")[, 4],F_Fledgling_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_fledgling_survival")[, 4],F_Adult_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_adult_survival")[, 4],M_Nestling_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_nestling_survival")[, 4],M_Groundling_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_groundling_survival")[, 4],M_Fledgling_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_fledgling_survival")[, 4],M_Adult_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_adult_survival")[, 4],Egg_survival =filter(survival_rates_boot_summary, species == species_name & vital_rate =="NA_egg_survival")[, 4],# Define h (harem size, h = 1 is monogamy) and k (clutch size)h = h,k = k,# Define primary sex ratioHSR = HSR) }
make_mprime_matrix()
This function makes a prime-matrix (i.e., a matrix halfway between the treatment matrix and the unbiased matrix) from a set of vital rate summaries
arguments:
survival_rates_boot_summary: the $vital_rate_ests_boot element in the tidied _hazard_rate_boot_tidy object from above (a summary of the bootstrapped vital rates
species: species name (either “BC” or “WBC”)
h: the mating system index specified in the parameter_distributions
k: the clutch size specified in the parameter_distributions
HSR: the hatching sex ratio specified in the parameter_distributions
sex: the sex that you wish to specify as the baseline vital rates for the LTRE
Code
make_mprime_matrix <-function(survival_rates_boot_summary, species_name, h, k, HSR, sex){if(sex =="male"){list(F_Nestling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_nestling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_nestling_survival")[, 4]) /2,F_Groundling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_groundling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_groundling_survival")[, 4]) /2,F_Fledgling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_fledgling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_fledgling_survival")[, 4]) /2,F_Adult_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_adult_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_adult_survival")[, 4]) /2,M_Nestling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_nestling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_nestling_survival")[, 4]) /2,M_Groundling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_groundling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_groundling_survival")[, 4]) /2,M_Fledgling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_fledgling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_fledgling_survival")[, 4]) /2,M_Adult_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_adult_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_adult_survival")[, 4]) /2,Egg_survival =filter(survival_rates_boot_summary, species == species_name)[15, 4],# Define h (harem size, h = 1 is monogamy) and k (clutch size)h = (h +1) /2,k = k,# Define primary sex ratioHSR = (HSR +0.5) /2) }else{list(F_Nestling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_nestling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_nestling_survival")[, 4]) /2,F_Groundling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_groundling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_groundling_survival")[, 4]) /2,F_Fledgling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_fledgling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_fledgling_survival")[, 4]) /2,F_Adult_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_adult_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_adult_survival")[, 4]) /2,M_Nestling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_nestling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_nestling_survival")[, 4]) /2,M_Groundling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_groundling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_groundling_survival")[, 4]) /2,M_Fledgling_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_fledgling_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_fledgling_survival")[, 4]) /2,M_Adult_survival = (filter(survival_rates_boot_summary, species == species_name & vital_rate =="Female_adult_survival")[, 4] +filter(survival_rates_boot_summary, species == species_name & vital_rate =="Male_adult_survival")[, 4]) /2,Egg_survival =filter(survival_rates_boot_summary, species == species_name)[15, 4],# Define h (harem size, h = 1 is monogamy) and k (clutch size)h = (h +1) /2,k = k,# Define primary sex ratioHSR = (HSR +0.5) /2) } }
LTRE_analysis()
LTRE_analysis() compares the relative difference in ASR response of an unbiased two-sex matrix vs. the observed sex-specific matrix
arguments:
Mprime_sens: the object produced by the sensitivity_analysis() applied to the Mprime matrix
matrix_str: the object of the matrix_structure expression
vital_rates: the object produced by the make_treat_matrix()
species_name: species name (either “BC” or “WBC”)
sex: the sex that you wish to specify as the baseline vital rates for the LTRE
Here we conduct the Life Table Response Experiment to decompose the difference in response between two or more “treatments” by weighting the difference in parameter values by the parameter’s contribution to the response (i.e., its sensitivity), and summing over all parameters. Our LTRE compares the observed scenario (M), to a null scenario with no sex differences (M0), whereby all male survival rates were set equal to the female rates (M0♀), the hatching sex ratio was unbiased (i.e., ρ = 0.5), and mating system was monogamous (i.e., h = 1). Thus, our LTRE identifies the drivers of ASR bias by decomposing the absolute parameter contributions to the difference between the ASR predicted by our model and an unbiased ASR.
Code
BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <-as.factor(BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <-as.factor(WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)# summarize the vital ratessurvival_rates_boot_summary <-bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>%mutate(vital_rate =paste(sex, stage, rate, sep ="_")) %>% Rmisc::summarySE(.,measurevar ="value",groupvars =c("vital_rate", "species"),conf.interval =0.95) %>%arrange(species, vital_rate)BC_VR_treat <-make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,h =1/pull(filter(parameter_distributions, species =="BC"& trait =="mating_system"), mean),HSR =pull(filter(parameter_distributions, species =="BC"& trait =="hatching_sex_ratio"), mean),k =pull(filter(parameter_distributions, species =="BC"& trait =="clutch_size"), mean),species_name ="BC")BC_VR_mprime_male <-make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,h =1/pull(filter(parameter_distributions, species =="BC"& trait =="mating_system"), mean),HSR =pull(filter(parameter_distributions, species =="BC"& trait =="hatching_sex_ratio"), mean),k =pull(filter(parameter_distributions, species =="BC"& trait =="clutch_size"), mean),species_name ="BC",sex ="male")BC_VR_mprime_female <-make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,h =1/pull(filter(parameter_distributions, species =="BC"& trait =="mating_system"), mean),HSR =pull(filter(parameter_distributions, species =="BC"& trait =="hatching_sex_ratio"), mean),k =pull(filter(parameter_distributions, species =="BC"& trait =="clutch_size"), mean),species_name ="BC",sex ="female")WBC_VR_treat <-make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,h =1/pull(filter(parameter_distributions, species =="WBC"& trait =="mating_system"), mean),HSR =pull(filter(parameter_distributions, species =="WBC"& trait =="hatching_sex_ratio"), mean),k =pull(filter(parameter_distributions, species =="WBC"& trait =="clutch_size"), mean),species_name ="WBC")WBC_VR_mprime_male <-make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,h =1/pull(filter(parameter_distributions, species =="WBC"& trait =="mating_system"), mean),HSR =pull(filter(parameter_distributions, species =="WBC"& trait =="hatching_sex_ratio"), mean),k =pull(filter(parameter_distributions, species =="WBC"& trait =="clutch_size"), mean),species_name ="WBC",sex ="male")WBC_VR_mprime_female <-make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,h =1/pull(filter(parameter_distributions, species =="WBC"& trait =="mating_system"), mean),HSR =pull(filter(parameter_distributions, species =="WBC"& trait =="hatching_sex_ratio"), mean),k =pull(filter(parameter_distributions, species =="WBC"& trait =="clutch_size"), mean),species_name ="WBC",sex ="female")BC_treatment_matrix <-coucal_matrix(BC_VR_treat)BC_M_prime_matrix_male <-coucal_matrix(BC_VR_mprime_male)BC_M_prime_matrix_female <-coucal_matrix(BC_VR_mprime_female)WBC_treatment_matrix <-coucal_matrix(WBC_VR_treat)WBC_M_prime_matrix_male <-coucal_matrix(WBC_VR_mprime_male)WBC_M_prime_matrix_female <-coucal_matrix(WBC_VR_mprime_female)BC_treatment_ASR_analysis <-matrix_ASR(M = BC_treatment_matrix, h = BC_VR_treat$h, k = BC_VR_treat$k,HSR = BC_VR_treat$HSR,iterations =100, num_boot =1, iter_add =1,species ="BC")BC_ASR_treat <- BC_treatment_ASR_analysis$ASRWBC_treatment_ASR_analysis <-matrix_ASR(M = WBC_treatment_matrix, h = WBC_VR_treat$h, k = WBC_VR_treat$k,HSR = WBC_VR_treat$HSR,iterations =100, num_boot =1, iter_add =1,species ="BC")WBC_ASR_treat <- WBC_treatment_ASR_analysis$ASRBC_M_prime_ASR_analysis_male <-matrix_ASR(M = BC_M_prime_matrix_male, h = BC_VR_mprime_male$h, k = BC_VR_mprime_male$k,HSR = BC_VR_mprime_male$HSR,iterations =100,num_boot =1,iter_add =1,species ="BC")BC_ASR_mprime_male <- BC_M_prime_ASR_analysis_male$ASRBC_M_prime_ASR_analysis_female <-matrix_ASR(M = BC_M_prime_matrix_female, h = BC_VR_mprime_female$h, k = BC_VR_mprime_female$k,HSR = BC_VR_mprime_female$HSR,iterations =100,num_boot =1,iter_add =1,species ="BC")BC_ASR_mprime_female <- BC_M_prime_ASR_analysis_female$ASRWBC_M_prime_ASR_analysis_male <-matrix_ASR(M = WBC_M_prime_matrix_male, h = WBC_VR_mprime_male$h, k = WBC_VR_mprime_male$k,HSR = WBC_VR_mprime_male$HSR,iterations =100,num_boot =1,iter_add =1,species ="BC")WBC_ASR_mprime_male <- WBC_M_prime_ASR_analysis_male$ASRWBC_M_prime_ASR_analysis_female <-matrix_ASR(M = WBC_M_prime_matrix_female, h = WBC_VR_mprime_female$h, k = WBC_VR_mprime_female$k,HSR = WBC_VR_mprime_female$HSR,iterations =100,num_boot =1,iter_add =1,species ="BC")WBC_ASR_mprime_female <- WBC_M_prime_ASR_analysis_female$ASRBC_lambda_treat <- BC_treatment_ASR_analysis$lambdaBC_lambda_mprime_male <- BC_M_prime_ASR_analysis_male$lambdaBC_lambda_mprime_female <- BC_M_prime_ASR_analysis_female$lambdaWBC_lambda_treat <- WBC_treatment_ASR_analysis$lambdaWBC_lambda_mprime_male <- WBC_M_prime_ASR_analysis_male$lambdaWBC_lambda_mprime_female <- WBC_M_prime_ASR_analysis_female$lambdaBC_treat_sensitivity_analysis <-sensitivity_analysis(vital_rate_summary = BC_VR_treat, matrix_str = matrix_structure, h = BC_VR_treat$h, k = BC_VR_treat$k, HSR = BC_VR_treat$HSR, niter =100, ASR = BC_ASR_treat,lambda = BC_lambda_treat)BC_Mprime_sensitivity_analysis_male <-sensitivity_analysis(vital_rate_summary = BC_VR_mprime_male, matrix_str = matrix_structure, h = BC_VR_mprime_male$h, k = BC_VR_mprime_male$k, HSR = BC_VR_mprime_male$HSR, niter =100, ASR = BC_ASR_mprime_male,lambda = BC_lambda_mprime_male)BC_Mprime_sensitivity_analysis_female <-sensitivity_analysis(vital_rate_summary = BC_VR_mprime_female, matrix_str = matrix_structure, h = BC_VR_mprime_female$h, k = BC_VR_mprime_female$k, HSR = BC_VR_mprime_female$HSR, niter =100, ASR = BC_ASR_mprime_female,lambda = BC_lambda_mprime_female) WBC_treat_sensitivity_analysis <-sensitivity_analysis(vital_rate_summary = WBC_VR_treat, matrix_str = matrix_structure, h = WBC_VR_treat$h, k = WBC_VR_treat$k, HSR = WBC_VR_treat$HSR, niter =100, ASR = WBC_ASR_treat,lambda = WBC_lambda_treat)WBC_Mprime_sensitivity_analysis_male <-sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_male, matrix_str = matrix_structure, h = WBC_VR_mprime_male$h, k = WBC_VR_mprime_male$k, HSR = WBC_VR_mprime_male$HSR, niter =100, ASR = WBC_ASR_mprime_male,lambda = WBC_lambda_mprime_male) WBC_Mprime_sensitivity_analysis_female <-sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_female, matrix_str = matrix_structure, h = WBC_VR_mprime_female$h, k = WBC_VR_mprime_female$k, HSR = WBC_VR_mprime_female$HSR, niter =100, ASR = WBC_ASR_mprime_female,lambda = WBC_lambda_mprime_female)
# Determine how much larger the contribution of each vital rates is compared to juvenile survival# groundling vs nestling:LTRE_coucal_ASR_prop <-filter(LTRE_coucal_ASR) %>% dplyr::select(-model) %>%mutate(parameter2 =as.character(parameter)) %>% dplyr::rename(parameter1 = parameter)prop_contribution_table <- LTRE_coucal_ASR_prop %>% tidyr::expand(., parameter1 = parameter1, parameter2 = parameter1) %>%left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter1, contribution, sex, species), by =c("parameter1")) %>% dplyr::rename(contribution1 = contribution) %>%left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter2, contribution, sex, species), by =c("parameter2", "species", "sex")) %>% dplyr::rename(contribution2 = contribution) %>%mutate(prop_diff =abs(contribution1) /abs(contribution2),prop_diff2 =ifelse(contribution1 <0& contribution2 <0, -(contribution1/contribution2),ifelse(contribution1 <0& contribution2 >0, contribution1/contribution2,ifelse(contribution1 >0& contribution2 <0, -(contribution1/contribution2), contribution1/contribution2))),parameter1 =factor(parameter1, levels =c("Hatching sex ratio","Nestling survival","Groundling survival","Fledgling survival","Adult survival","Mating system")),parameter2 =factor(parameter2, levels =c("Hatching sex ratio","Nestling survival","Groundling survival","Fledgling survival","Adult survival","Mating system")))LTRE_heatmap <- prop_contribution_table %>%filter(parameter1 %in%c("Mating system","Adult survival","Fledgling survival","Groundling survival","Nestling survival","Hatching sex ratio") & parameter2 %in%c("Mating system","Adult survival","Fledgling survival","Groundling survival","Nestling survival","Hatching sex ratio")) %>%ggplot(., aes(parameter1, parameter2, fill = prop_diff2)) +geom_tile() +facet_grid(sex ~ species,labeller =labeller(.cols = species_names, .rows = analysis_names)) +theme(axis.text.x =element_text(size =8, angle =45, hjust =1),axis.text.y =element_text(size =8),axis.title.x =element_text(size =12),axis.title.y =element_text(size =12),legend.position ="bottom",legend.box ="horizontal",legend.text =element_text(size =7),legend.title =element_text(size =7),axis.ticks =element_blank(),panel.border =element_rect(color ="grey40", fill =NA),strip.background =element_blank(),strip.text =element_text(size =12, colour ="black", face ="italic")) +scale_x_discrete(labels =c("Hatching sex ratio"=expression(italic("\u03C1")),"Nestling survival"=expression(italic(S["n"])),"Groundling survival"=expression(italic(S["g"])),"Fledgling survival"=expression(italic(S["f"])),"Adult survival"=expression(italic(phi["a"])),"Mating system"=expression(italic("h")))) +scale_y_discrete(labels =c("Hatching sex ratio"=expression(italic("\u03C1")),"Nestling survival"=expression(italic(S["n"])),"Groundling survival"=expression(italic(S["g"])),"Fledgling survival"=expression(italic(S["f"])),"Adult survival"=expression(italic(phi["a"])),"Mating system"=expression(italic("h")))) +scale_fill_gradient2(low = sex_pal2[1],mid ="white",high = sex_pal2[2], name ="Proportional contribution of A compared to B") +#,xlab("Parameter A") +ylab("Parameter B")LTRE_heatmap
Supplementary Analysis 1: sex-specific juvenile movement
Here we assess individual heterogeneity in several movement traits based on radio telemetry data of tagged juveniles
Code
coucal_wp <-read.csv("data/raw/juvies_2013_2020_waypoints.csv",stringsAsFactors =FALSE, header =TRUE) %>%# clean up a few inconsistencies in the datamutate(timestamp =ifelse(ring_ID =="GN 76660"& date_dec =="20140502", "2014-05-02 10:09:00.000", timestamp),location.long =ifelse(location.long ==34.68357, 34.08357, location.long)) %>%filter(ring_ID !="GN 76726") %>%filter(ring_ID !="GN 76728") %>%filter(ring_ID !="GN 76732") %>%filter(ring_ID !="GN 76733") %>%filter(ring_ID !="GN 76777") %>%# both male and female??filter(movebank_event.id !="16897537877") %>%filter(movebank_event.id !="16897535798") %>%filter(movebank_event.id !="16897537758") %>%filter(movebank_event.id !="16897538205") %>%# specify the timezone of the time stampmutate(timestamp =as.POSIXct(timestamp,tz ="Africa/Nairobi"),# make the serial number a 4 digit stringserial_number =as.character(str_pad(serial_number, 4, pad ="0"))) %>%# make a julian date time variable for plottingmutate(julian_date =as.numeric(format(timestamp, "%j"))) %>%# remove observations older than 70 days of agefilter(age_days <71) %>%# consolidate columns of interest dplyr::select(species, year, nest.nr, ring_ID, sex, timestamp, julian_date, age_days, location.long, location.lat, species, site, location, comment)# specify coordinate columnscoordinates(coucal_wp) <-c("location.long","location.lat")# define spatial projection as WGS84proj4string(coucal_wp) <-CRS("+init=epsg:4326")
Explore the spatial extent of the data
Code
# specify mapview options for viewingmapviewOptions(basemaps =c("Esri.WorldImagery"),layers.control.pos ="topright",legend.pos ="bottomleft")# check out the spatial datamapview(coucal_wp[coucal_wp$ring_ID %in%"GN 76751",], zcol ="julian_date",layer.name ="WBC radio tagging",layers.control.pos ="topright")
Set-up the spatial data and create a trajectory format for analysis
Code
# make a spatial trajectory object of each animal's encounter historytr_coucal_wp <-as.ltraj(xy = sp::coordinates(coucal_wp),date = coucal_wp$timestamp,id = coucal_wp$ring_ID)# convert the trajectory object back to a dataframe, rename columns, and# calculate the cumulative distance traveled over the observationtr_coucal_wp <-ld(tr_coucal_wp) %>% dplyr::rename(distance = dist,dispersion = R2n,cardinal_angle = abs.angle,relative_angle = rel.angle,ring_ID = id) %>%group_by(ring_ID) %>%mutate(cum_distance =cumsum(distance))# extract the temporal summary info for each bird (min age, min date, max age) # and merge to the trajectory dataframecoucal_wp_clean <-as.data.frame(coucal_wp) %>%group_by(species, ring_ID, sex, nest.nr) %>% dplyr::summarise(min_age =min(age_days),min_date =min(timestamp),max_age =max(age_days)) %>%arrange(desc(max_age)) %>%left_join(., tr_coucal_wp, by ="ring_ID") %>%# calcuate the differences in days between the observation and minimum date# to get time since first observation, convert from seconds to days, then add# to the minimum age value to extract the age at a given observationmutate(date_diff = date - min_date,year =year(date)) %>%mutate(date_diff =round(((as.numeric(date_diff)/60)/60)/24)) %>%mutate(age = min_age + date_diff) %>%ungroup() %>%# remove all NA rowsfilter(!is.na(cum_distance)) %>%# calculate the temporal lag between observationsgroup_by(ring_ID) %>%mutate(obs_diff = date_diff -lag(date_diff),origin_x = x[which.min(date)],origin_y = y[which.min(date)])# for each bird calculate the daily distance from origincoucal_wp_clean$distance_origin <-sapply(1:nrow(coucal_wp_clean), function(i)spDistsN1(as.matrix(coucal_wp_clean[i, c("x", "y")]), as.matrix(coucal_wp_clean[i, c("origin_x", "origin_y")]), longlat =TRUE))# transform the distance metrics to meterscoucal_wp_clean <- coucal_wp_clean %>%mutate(distance = distance *100000,distance_origin = distance_origin *1000,cum_distance = cum_distance *100000)# check the distribution of distanceshist(log(coucal_wp_clean$distance))
# extract individuals that mainly had observations spaced 2 days aparttwo_day_rings <- coucal_wp_clean %>% dplyr::select(species, ring_ID, sex, date_diff, obs_diff, age) %>%group_by(ring_ID) %>% dplyr::summarise(mean_obs_diff =mean(obs_diff, na.rm =TRUE),median_obs_diff =median(obs_diff, na.rm =TRUE)) %>%arrange(desc(median_obs_diff)) %>%filter(median_obs_diff <=2)# subset to only include individuals with mainly observations spaced 2 days apartcoucal_wp_clean <- coucal_wp_clean %>%filter(ring_ID %in% two_day_rings$ring_ID) %>%# remove questionable observation 18 km filter(distance <18000)
Model 1: daily distance moved
Does daily distance moved increase with age since leaving the nest and does this relationship differ between males and females?
Code
#### Modeling sex-specific movements ##### M1 Black Coucals: strong effect of age, no interaction with sexmod_BC_dist <-lmer(log(distance) ~poly(age, 2) * sex + (1|ring_ID),data =filter(coucal_wp_clean, species =="BC"))# effect-size plotscoefplot::coefplot(mod_BC_dist)
Code
plot(allEffects(mod_BC_dist))
Code
# extract fitted valuesmod_BC_dist_fits <-as.data.frame(effect(term ="poly(age, 2) * sex", mod = mod_BC_dist, xlevels =list(age =seq(min(coucal_wp_clean[coucal_wp_clean$species =="BC", "age"]),max(coucal_wp_clean[coucal_wp_clean$species =="BC", "age"]), 1),sex =c("Female", "Male"))))mod_BC_dist_fits$species <-"BC"# M1 White-browed Coucals: strong effect of age, no interaction with sexmod_WBC_dist <-lmer(log(distance) ~poly(age, 2) * sex + (1|ring_ID), data =filter(coucal_wp_clean, species =="WBC"& distance >0))# effect-size plotscoefplot::coefplot(mod_WBC_dist)
# M1 plotggplot() + luke_theme +theme(legend.background =element_rect(fill="transparent")) +geom_line(data = coucal_wp_clean, aes(x = age, y =log(distance), group = ring_ID, color = sex),alpha =0.2) +geom_line(data = mod_dist_fits, aes(x = age, y = fit, color = sex),lwd =0.5) +geom_ribbon(data = mod_dist_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),lwd =1, alpha =0.25) +facet_grid(species ~ ., labeller =labeller(species = species.labs)) +ylab(expression(paste("Distance moved between observations (log10)"%+-%"95% CI", sep =""))) +xlab("Age since leaving nest") +scale_color_manual(values = sex_pal2,name ="Sex",breaks =c("Female", "Male"),labels =c("Female", "Male")) +scale_fill_manual(values = sex_pal2,name ="Sex",breaks =c("Female", "Male"),labels =c("Female", "Male"))
Code
# # first set up plotting parameters# # define the plotting theme to be used in subsequent ggplots# luke_theme <- # theme_bw() +# theme(# text = element_text(family = "Verdana"),# legend.title = element_text(size = 16),# legend.text = element_text(size = 12),# axis.title.x = element_text(size = 12),# axis.text.x = element_text(size = 12), # axis.title.y = element_text(size = 16),# axis.text.y = element_text(size = 12),# strip.text = element_text(size = 12),# panel.grid.major = element_blank(),# panel.grid.minor = element_blank(),# axis.ticks = element_line(size = 0.5, colour = "grey40"),# axis.ticks.length = unit(0.2, "cm"),# panel.border = element_rect(linetype = "solid", colour = "grey"),# legend.position = c(0.1, 0.9)# )# # # set plotting color palettes# sex_pal2 <- RColorBrewer::brewer.pal(8, "Dark2")[c(2,1)]# plot_palette_brood <- RColorBrewer::brewer.pal(7, "Blues")[c(3:7)]# # # specify the facet labels for each species
Model 2: cumulative distance moved
Does cumulative distance moved increase with age since leaving the nest and does this relationship differ between males and females?
Code
# M2 Black Coucals: strong effect of age, no interaction with sexmod_BC_cum_dist <-lmer(log(cum_distance) ~poly(age, 2) * sex + (1|ring_ID),data =filter(coucal_wp_clean, species =="BC"))# effect-size plotscoefplot::coefplot(mod_BC_cum_dist)
Code
plot(allEffects(mod_BC_cum_dist))
Code
# extract fitted valuesmod_BC_cum_dist_fits <-as.data.frame(effect(term ="poly(age, 2) * sex", mod = mod_BC_cum_dist, xlevels =list(age =seq(min(coucal_wp_clean[coucal_wp_clean$species =="BC", "age"]),max(coucal_wp_clean[coucal_wp_clean$species =="BC", "age"]), 1),sex =c("Female", "Male"))))mod_BC_cum_dist_fits$species <-"BC"# M2 White-browed Coucals: strong effect of age, no interaction with sexmod_WBC_cum_dist <-lmer(log(cum_distance) ~poly(age, 2) * sex + (1|ring_ID), data =filter(coucal_wp_clean, species =="WBC"& distance >0))# effect-size plotscoefplot::coefplot(mod_WBC_cum_dist)
Do individuals move further from the nest as they age and does this relationship differ between males and females?
Code
# M3 Black Coucals: strong effect of age, no interaction with sexmod_BC_disp <-lmer(log(distance_origin) ~poly(age, 2) * sex + (1|ring_ID),data =filter(coucal_wp_clean, species =="BC"& distance_origin >0))# effect-size plotscoefplot::coefplot(mod_BC_disp)
Code
plot(allEffects(mod_BC_disp))
Code
# extract fitted valuesmod_BC_disp_fits <-as.data.frame(effect(term ="poly(age, 2) * sex", mod = mod_BC_disp, xlevels =list(age =seq(min(coucal_wp_clean[coucal_wp_clean$species =="BC", "age"]),max(coucal_wp_clean[coucal_wp_clean$species =="BC", "age"]), 1),sex =c("Female", "Male"))))mod_BC_disp_fits$species <-"BC"# M3 White-browed Coucals: strong effect of age, no interaction with sexmod_WBC_disp <-lmer(log(distance_origin) ~poly(age, 2) * sex + (1|ring_ID), data =filter(coucal_wp_clean, species =="WBC"& distance_origin >0))# effect-size plotscoefplot::coefplot(mod_WBC_disp)
Here we run a weekly mark-recapture survival analysis on the weekly encounters (and dead recoveries) of adults throughout the season. The first step is to wrangle the capture history into the format needed for a Burnham-joint-live-dead model: each encounter event requires two characters (00 = not encountered, 01 = encountered alive, 11 = encountered dead)
Code
status_dat_ad <-# read raw dataread.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv", header =TRUE, stringsAsFactors =FALSE, na.strings =c("", " ", "NA")) %>%# rename ring_ID column dplyr::rename(ring_ID = Alu,ageC2 = days_from_ringing_until_last_observation) %>%# specify date strings properlymutate(month_ringed =str_sub(date_dec_initially_ringed, start =5, end =6),day_ringed =str_sub(date_dec_initially_ringed, start =7, end =8),year_ringed =str_sub(date_dec_initially_ringed, start =1, end =4),month_last_ob =str_sub(date_dec_last_observed, start =5, end =6),day_last_ob =str_sub(date_dec_last_observed, start =7, end =8),year_last_ob =str_sub(date_dec_initially_ringed, start =1, end =4)) %>%mutate(date_ringed =as.Date(paste(year_ringed, month_ringed, day_ringed, sep ="-"), format ="%Y-%m-%d"),date_last_ob =as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep ="-"), format ="%Y-%m-%d")) %>%# select variables of interest dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>%# remove all white space from datamutate(across(everything(), ~str_trim(.x))) %>%mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%# specify empty data as NAmutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%# exclude all individuals that died in the nest# filter(Fledged_status == "yes") %>% # classify columnsmutate(sex =as.factor(sex),ageC2 =as.numeric(ageC2)) %>%# remove rows with missing sex, age, and status infofilter(!is.na(sex) &!is.na(ageC2) &!is.na(statusC)) %>%# make a unique id for each individualmutate(# create the age of entry into the data (all at age 15)entry =0,# specify the age of death or censoringexit = ageC2,# make the event numeric and specify if # the individual died (1) or was censored (0)event =as.numeric(statusC),sex =ifelse(str_detect(sex, pattern ="[Ff]emale"), "F", ifelse(str_detect(sex, pattern ="[Mm]ale"), "M", sex))) %>%# consolidate to variables of interest dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob)status_dat_ad_week <- status_dat_ad %>% dplyr::select(species, ring_ID, sex, date_ringed, date_last_ob) %>%pivot_longer(-c(species, ring_ID, sex), names_to ="status", values_to ="date") %>%mutate(year_numeric =as.numeric(lubridate::year(date)),week_numeric =as.numeric(strftime(date, format ="%V"))) %>% dplyr::select(species, ring_ID, year_numeric, sex, week_numeric)# import raw csv data into RBC_detect_dat_A <-read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na ="NA", col_types ="text") %>% dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>%mutate(month =str_sub(date_dec, start =5, end =6),day =str_sub(date_dec, start =7, end =8)) %>%mutate(date =as.Date(paste(year, month, day, sep ="-"), format ="%Y-%m-%d")) %>%mutate(across(c("sex", "site", "age_status"), tolower)) %>%mutate(age_status =ifelse(age_status =="adult", "A", ifelse(age_status =="juvenile", "J", "XXX")),ring_ID =str_replace_all(string = ring_ID, fixed(" "), "")) %>%mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>%mutate(sex =ifelse(str_detect(sex, pattern ="[Ff]emale"), "F",ifelse(str_detect(sex, pattern ="[Mm]ale"), "M", sex)),month_year =format(date, "%Y-%m"),month_numeric =as.numeric(month),year_numeric =as.numeric(year),week_numeric =as.numeric(strftime(date, format ="%V"))) %>%filter(species =="BC") %>%filter(age_status =="A") %>% dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%bind_rows(., filter(status_dat_ad_week, species =="BC")) %>%filter(!is.na(week_numeric)) %>%distinct() %>%mutate(week_std =round(scale_by(week_numeric ~ year_numeric, ., scale =0), digits =0)) %>%mutate(week_std = week_std +abs(min(week_std, na.rm =TRUE))) %>% dplyr::rename(year = year_numeric)WBC_detect_dat_A <-read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na ="NA", col_types ="text") %>% dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>%mutate(month =str_sub(date_dec, start =5, end =6),day =str_sub(date_dec, start =7, end =8)) %>%mutate(date =as.Date(paste(year, month, day, sep ="-"), format ="%Y-%m-%d")) %>%mutate(across(c("sex", "site", "age_status"), tolower)) %>%mutate(age_status =ifelse(age_status =="adult", "A", ifelse(age_status =="juvenile", "J", "XXX")),ring_ID =str_replace_all(string = ring_ID, fixed(" "), "")) %>%mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>%mutate(sex =ifelse(str_detect(sex, pattern ="[Ff]emale"), "F",ifelse(str_detect(sex, pattern ="[Mm]ale"), "M", sex)),month_year =format(date, "%Y-%m"),month_numeric =as.numeric(month),year_numeric =as.numeric(year),week_numeric =as.numeric(strftime(date, format ="%V"))) %>%filter(species =="WBC") %>%filter(age_status =="A") %>% dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%bind_rows(., filter(status_dat_ad_week, species =="WBC")) %>%filter(!is.na(week_numeric)) %>%distinct() %>%mutate(week_std =round(scale_by(week_numeric ~ year_numeric, ., scale =0), digits =0)) %>%mutate(week_std = week_std +abs(min(week_std, na.rm =TRUE))) %>% dplyr::rename(year = year_numeric)# use the BaSTA function "CensusToCaptHist()" function to convert long format# encounter histories of each individual, to wide format with 1's and 0's for # each week of encounter for each yearBC_weekly_A_ch <-CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,d = BC_detect_dat_A$week_std,dformat ="%W", timeInt ="%W") %>% dplyr::rename(ring_ID = ID) %>%mutate(ID =rownames(.)) %>%left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by ="ring_ID") %>%distinct()WBC_weekly_A_ch <-CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,d = WBC_detect_dat_A$week_std,dformat ="%W", timeInt ="%W") %>% dplyr::rename(ring_ID = ID) %>%mutate(ID =rownames(.)) %>%left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by ="ring_ID") %>%distinct()# create a two-character string for each encounter and clean the outputBC_weekly_A_ch_Burnham <-sapply(BC_weekly_A_ch[2:26], function(x) paste(x, "0", sep ="")) %>%cbind(BC_weekly_A_ch[c(1,45:length(BC_weekly_A_ch))]) %>%mutate(across(everything(), ~str_replace(string = .x, pattern ="NA0", "00"))) %>%mutate(across(everything(), ~str_trim(.x))) %>%mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%distinct() %>%mutate(sex =as.factor(sex))WBC_weekly_A_ch_Burnham <-sapply(WBC_weekly_A_ch[2:26], function(x) paste(x, "0", sep ="")) %>%cbind(WBC_weekly_A_ch[c(1,54:length(WBC_weekly_A_ch))]) %>%mutate(across(everything(), ~str_replace(string = .x, pattern ="NA0", "00"))) %>%mutate(across(everything(), ~str_trim(.x))) %>%mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%distinct() %>%mutate(sex =as.factor(sex))# Import status datastatus_dat_all_ad <-# read raw dataread.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv", header =TRUE, stringsAsFactors =FALSE, na.strings =c("", " ", "NA")) %>%# rename ring_ID column dplyr::rename(ring_ID = Alu,ageC2 = days_from_ringing_until_last_observation) %>%# specify date strings properlymutate(month_ringed =str_sub(date_dec_initially_ringed, start =5, end =6),day_ringed =str_sub(date_dec_initially_ringed, start =7, end =8),year_ringed =str_sub(date_dec_initially_ringed, start =1, end =4),month_last_ob =str_sub(date_dec_last_observed, start =5, end =6),day_last_ob =str_sub(date_dec_last_observed, start =7, end =8),year_last_ob =str_sub(date_dec_initially_ringed, start =1, end =4)) %>%mutate(date_ringed =as.Date(paste(year_ringed, month_ringed, day_ringed, sep ="-"), format ="%Y-%m-%d"),date_last_ob =as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep ="-"), format ="%Y-%m-%d")) %>%# select variables of interest dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>%# remove all white space from datamutate(across(everything(), ~str_trim(.x))) %>%mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%# specify empty data as NAmutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%# classify columnsmutate(sex =as.factor(sex),ageC2 =as.numeric(ageC2)) %>%# remove rows with missing sex, age, and status infofilter(!is.na(sex) &!is.na(ageC2) &!is.na(statusC)) %>%# make a unique id for each individualmutate(# create the age of entry into the data (all at age 15)entry =0,# specify the age of death or censoringexit = ageC2,# make the event numeric and specify if # the individual died (1) or was censored (0)event =as.numeric(statusC),sex =ifelse(str_detect(sex, pattern ="[Ff]emale"), "F", ifelse(str_detect(sex, pattern ="[Mm]ale"), "M", sex))) %>%# consolidate to variables of interest dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob) %>% dplyr::select(species, ring_ID, year, sex, event)# join detection history with status dataBC_detect_status_join <-filter(status_dat_all_ad, species =="BC") %>%left_join(BC_weekly_A_ch_Burnham, ., by =c("ring_ID", "year", "sex")) %>% dplyr::select(-species) %>%mutate(event =ifelse(is.na(event), 0, event))WBC_detect_status_join <-filter(status_dat_all_ad, species =="WBC") %>%left_join(WBC_weekly_A_ch_Burnham, ., by =c("ring_ID", "year", "sex")) %>% dplyr::select(-species) %>%mutate(event =ifelse(is.na(event), 0, event))# determine the last and first detection for each individual max_age_index_BC <-apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x =="10")) %>%lapply(., function(x) x[which.max(x)]) %>%unlist(.)min_age_index_BC <-apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x =="10")) %>%lapply(., function(x) x[which.min(x)]) %>%unlist(.)max_age_index_WBC <-apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x =="10")) %>%lapply(., function(x) x[which.max(x)]) %>%unlist(.)min_age_index_WBC <-apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x =="10")) %>%lapply(., function(x) x[which.min(x)]) %>%unlist(.)# put "11" at the last detection if the status was a "1" (i.e., dead)for(i in1:nrow(BC_detect_status_join)){if(BC_detect_status_join$event[i] =="1"){ BC_detect_status_join[i, max_age_index_BC[i]] <-"11" }}for(i in1:nrow(WBC_detect_status_join)){if(WBC_detect_status_join$event[i] =="1"){ WBC_detect_status_join[i, max_age_index_WBC[i]] <-"11" }}# bind the min and max age info to the capture history datadetect_status_BC <-bind_cols(BC_detect_status_join, max_age_index_BC) %>%bind_cols(., min_age_index_BC) %>% dplyr::rename(std_week_last_obs ='...31',std_week_first_obs ='...32')detect_status_WBC <-bind_cols(WBC_detect_status_join, max_age_index_WBC) %>%bind_cols(., min_age_index_WBC) %>% dplyr::rename(std_week_last_obs ='...31',std_week_first_obs ='...32')# consolidate capture histories for each sexBlack_Coucal_adult_weekly_Burnham_ch <-apply(detect_status_BC[, c(1:25)], 1, paste, collapse ="") %>%as.data.frame() %>% dplyr::rename(ch ='.') %>%mutate(ch =as.character(ch)) %>% dplyr::bind_cols(., detect_status_BC[, c(26:length(detect_status_BC))]) %>%filter(str_detect(ch, "11", negate =TRUE) |str_count(ch, "1") >2) %>% dplyr::select(ch, ring_ID, sex, year) %>%drop_na() %>%mutate(year =ifelse(year =="1006", "2006", year)) %>%mutate_at(vars(ring_ID, sex, year), as.factor) White_browed_Coucal_adult_weekly_Burnham_ch <-apply(detect_status_WBC[, c(1:25)], 1, paste, collapse ="") %>%as.data.frame() %>% dplyr::rename(ch ='.') %>%mutate(ch =as.character(ch)) %>% dplyr::bind_cols(., detect_status_WBC[, c(26:length(detect_status_WBC))]) %>%filter(str_detect(ch, "11", negate =TRUE) |str_count(ch, "1") >2) %>% dplyr::select(ch, ring_ID, sex, year) %>%drop_na() %>%mutate_at(vars(ring_ID, sex, year), as.factor)BC_adult_ch <- Black_Coucal_adult_weekly_Burnham_ch %>%mutate(sex =as.factor(sex),year =as.factor(year)) %>% dplyr::select(ch, sex, year)WBC_adult_ch <- White_browed_Coucal_adult_weekly_Burnham_ch %>%mutate(sex =as.factor(sex),year =as.factor(year)) %>% dplyr::select(ch, sex, year)
Run the survival analysis in RMark, extract the model predictions, and plot the results.
Code
# process the adult data as a "Burnham" analysisBC_coucal_adult.proc <- RMark::process.data(data = BC_adult_ch,model ="Burnham",groups =c("sex", "year"))WBC_coucal_adult.proc <- RMark::process.data(data = WBC_adult_ch,model ="Burnham",groups =c("sex", "year"))# make design matrixBC_coucal_adult.ddl <- RMark::make.design.data(BC_coucal_adult.proc)WBC_coucal_adult.ddl <- RMark::make.design.data(WBC_coucal_adult.proc)burnham_S_analysis_run_pFr_dot <-function(proc_data, design_data){# sex-specific model S.sex =list(formula =~sex)# sex- and week-specific model S.sex_Time =list(formula =~sex + Time)# sex- and week-specific model (quadratic) S.sex_Quad =list(formula =~sex + (I(Time) +I(Time^2)))# sex- and week-specific model (interaction) S.sex_x_Time =list(formula =~sex * Time)# sex- and week-specific model (quadratic with interaction) S.sex_x_Quad =list(formula =~sex * (I(Time) +I(Time^2)))# p (encounter probability):# constant model p.dot =list(formula =~1)# F (site fidelity probability):# constant model F.dot =list(formula =~1)# r (recovery probability)# constant model r.dot =list(formula =~1)# create model list for all a priori models above cml <- RMark::create.model.list("Burnham")# run the models in program MARK model.list <- RMark::mark.wrapper(cml, data = proc_data, ddl = design_data, delete =TRUE, wrap =FALSE, threads =1, brief =TRUE,silent =TRUE, output =FALSE)# output the model list and sotre the resultsreturn(model.list) }adult_S_analysis_out_BC <-burnham_S_analysis_run_pFr_dot(proc_data = BC_coucal_adult.proc,design_data = BC_coucal_adult.ddl)adult_S_analysis_out_WBC <-burnham_S_analysis_run_pFr_dot(proc_data = WBC_coucal_adult.proc,design_data = WBC_coucal_adult.ddl)BC_sex_x_quad_reals <- adult_S_analysis_out_BC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>%bind_cols(data.frame(str_split_fixed(rownames(.), " ", n =5)), .) %>%mutate(age =as.integer(unlist(str_extract_all(X4,"[0-9]+"))),sex =as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) =="F", "Female", "Male"))) %>%filter(X1 =="S") %>% dplyr::select(sex, age, estimate, se, lcl, ucl) %>%mutate(species ="Black Coucal")WBC_sex_x_quad_reals <- adult_S_analysis_out_WBC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>%bind_cols(data.frame(str_split_fixed(rownames(.), " ", n =5)), .) %>%mutate(age =as.integer(unlist(str_extract_all(X4,"[0-9]+"))),sex =as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) =="F", "Female", "Male"))) %>%filter(X1 =="S") %>% dplyr::select(sex, age, estimate, se, lcl, ucl) %>%mutate(species ="White-browed Coucal")BC_sex_reals <- adult_S_analysis_out_BC$S.sex.p.dot.r.dot.F.dot$results$real %>%bind_cols(data.frame(str_split_fixed(rownames(.), " ", n =5)), .) %>%mutate(age =as.integer(unlist(str_extract_all(X4,"[0-9]+"))),sex =as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) =="F", "Female", "Male"))) %>%filter(X1 =="S") %>% dplyr::select(sex, age, estimate, se, lcl, ucl) %>%mutate(species ="Black Coucal")WBC_sex_reals <- adult_S_analysis_out_WBC$S.sex.p.dot.r.dot.F.dot$results$real %>%bind_cols(data.frame(str_split_fixed(rownames(.), " ", n =5)), .) %>%mutate(age =as.integer(unlist(str_extract_all(X4,"[0-9]+"))),sex =as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) =="F", "Female", "Male"))) %>%filter(X1 =="S") %>% dplyr::select(sex, age, estimate, se, lcl, ucl) %>%mutate(species ="White-browed Coucal")time_varying_plot <-ggplot(data =bind_rows(BC_sex_x_quad_reals, WBC_sex_x_quad_reals)) +geom_line(aes(x = age, y = estimate, group = sex, color = sex), lwd =0.8) +geom_ribbon(aes(x = age, ymin = lcl, ymax = ucl, group = sex, fill = sex), alpha =0.5) +scale_color_manual(values = sex_pal2,labels =c("Female", "Male")) +scale_fill_manual(values = sex_pal2,labels =c("Female", "Male")) +facet_grid(. ~ species) +theme_bw() +theme(text =element_text(family ="Verdana"),axis.title =element_text(size =12),axis.text =element_text(size =10), strip.text =element_text(size =12),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_line(size =0.5, colour ="grey40"),axis.ticks.length =unit(0.2, "cm"),panel.border =element_rect(linetype ="solid", colour ="grey"),legend.position =c(0.8, 0.2),legend.title =element_blank()) +xlab("Week of breeding season (standardized across years)") +ylab("Adult weekly survival probability (± 95% CI)")constant_plot <- ggplot2::ggplot() +geom_errorbar(data =bind_rows(BC_sex_reals, WBC_sex_reals),aes(x = sex, ymax = ucl, ymin = lcl, color = sex),alpha =1, width =0.05, lwd =0.5) +geom_point(data =bind_rows(BC_sex_reals, WBC_sex_reals),aes(x = sex, y = estimate, fill = sex),lwd =3, shape =21) +scale_color_manual(values = sex_pal2,labels =c("Female", "Male")) +scale_fill_manual(values = sex_pal2,labels =c("Female", "Male")) +facet_grid(. ~ species) +theme_bw() +theme(text =element_text(family ="Verdana"),axis.title =element_blank(),axis.text =element_text(size =10), strip.text =element_text(size =12),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),axis.ticks =element_line(size =0.5, colour ="grey40"),axis.ticks.length =unit(0.2, "cm"),panel.border =element_rect(linetype ="solid", colour ="grey"),legend.position ="none") +ylab("Adult weekly survival probability (± 95% CI)") +scale_y_continuous(limits =c(0, 1))combo_plot_S <- time_varying_plot + constant_plot +plot_annotation(tag_levels ="A")combo_plot_S
Supplementary Analysis 3: Plumage-based sex-specific age structure
Here we assess annual variation in the sex-specific age structure of the Black Coucal population. Older individuals (i.e., 3+ years) have rufous flight feathers (i.e., primary coverts, primaries, secondary coverts, secondaries), whereas these flight feathers are barred for younger individuals
